Strongly nonlinear dynamics of electrolytes in large ac voltages 



o 
o 

< 

(N 

i 

q 

O 

•i-H 

>V 

Ph. 



O 

in 
rn 

o 

OS 

o 



Laurits H0jgaard Olesen, 1 Martin Z. Bazant, 2,3 and Henrik Bruus 1 

1 Department of Micro and Nanotechnology, Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark 
2 Departments of Chemical Engineering and Mathematics, 
Massachusetts Institute of Technology, Cambridge, MA 02139, USA 
3 UMR Gulliver ESPCI-CNRS 7083, 10 rue Vauquehn, F-75005 Pans, France 

(Dated: 24 August 2009) 

We study the response of a model micro-electrochemical cell to a large ac voltage of frequency 
comparable to the inverse cell relaxation time. To bring out the basic physics, we consider the 
simplest possible model of a symmetric binary electrolyte confined between parallel-plate blocking 
electrodes, ignoring any transverse instability or fluid flow. We analyze the resulting one-dimensional 
problem by matched asymptotic expansions in the limit of thin double layers and extend previous 
work into the strongly nonlinear regime, which is characterized by two novel features - significant 
salt depletion in the electrolyte near the electrodes and, at very large voltage, the breakdown of the 
quasi-equilibrium structure of the double layers. The former leads to the prediction of "ac capacitive 
desalination" , since there is a time-averaged transfer of salt from the bulk to the double layers, via 
oscillating diffusion layers. The latter is associated with transient diffusion limitation, which drives 
the formation and collapse of space-charge layers, even in the absence of any net Faradaic current 
through the cell. We also predict that steric effects of finite ion sizes (going beyond dilute solution 
theory) act to suppress the strongly nonlinear regime in the limit of concentrated electrolytes, ionic 
liquids and molten salts. Beyond the model problem, our reduced equations for thin double layers, 
based on uniformly valid matched asymptotic expansions, provide a useful mathematical framework 
to describe additional nonlinear responses to large ac voltages, such as Faradaic reactions, electro- 
osmotic instabilities, and induced-charge electrokinetic phenomena. 

PACS numbers: 



I. INTRODUCTION 

Time-dependent voltages are applied to electrolytes in 
many different fields, and theoretical models to interpret 
the results have been developed for over a century 
Current applications include energy storage in electro- 
chemical systems (e.g. supercapacitors high-rate 
batteries [1, [E0)> now control in microffuidics (e.g . ac 
electroosmotic EM, M El El El El El El, El tiM and 
electrothermal pjj, |2Q| flows), particle handling in col- 
loidal materials (e.g. dielectrophoresis [2ll. |22| . induced- 
charge electrophoresis [23, [24j, |25|, |26j), cellular and 
molecular manipulation in the biological systems (e.g . 
electroporation [2?], [H, l29ll. cell trappi ng [221 . [30l . l3lf . 
and biomolecular sensing [33. l33l . HI [35|, 

In many cases, periodic voltages are used to drive alter- 
nating current (ac) to eliminate any net linear response, 
such as direct current or electroosmotic flow. The most 
common application of ac forcing is in impedance spec- 
troscopy, long used to characterize electrochemical inter- 
faces [37j. The current response to a small sinusoidal 
voltage is fitted to an electrical circuit model, where the 
interface acts an impedance in series with a bulk re- 
sistance [H, HI 113] • The characteristic frequency for 
double-layer charging is then the inverse "i?C time" 
of the equivalent circuit [![. Circuit models are also 
used to describe electrochemical response in much more 
complicated situations, such as composite por ous elec- 
trodes [1, EJ, micro-electrode arrays d [ToL Il4. El, HH, 
and biological tissues [H HI] . 

Circuit models can be derived from underlying ion- 



transport equations by considering the joint limit of thin 
double layers (e = Xd/L <C 1, where A_d is the Debye- 
Hiickel screening length and L is the geometrical scale) 
and small voltages (V <C kT/e where V is the ampli- 
tude of the applied voltage and kT/e is the thermal volt- 
age) [l[. From a mathematical point of view, this can be 
done systematically starting from the Poisson-Nernst- 
Planck equations (PNP) by asymptotic boundary-layer 
analysis, which was introduced to electrochemistry in 
the 1960s to ju stify the thin-double-layer approxima- 
tion [HI, EI El H El- The joint asymptotic limit of 
thin double layers and large voltages, which is mathe- 
matically more challenging and physically more complex, 
has also been analyzed under conditions of steady direct 
current (dc). At sufficiently large steady dc currents, 
exceeding diffusion limitation, a variety of exotic effects 
arise, such as the expansion of the double layer into an 
extended, non-equilibrium space-charge layer [H El Ell 
and electro-hydrodynamic instability due to second-kind 
electroosmotic flows [H [5(| [5l|. Clearly, such effects 
cannot be captured by classical circuit models, but they 
continue to be used in dynamical situations, even with 
large voltages, for lack of a simple mathematical alterna- 
tive. 

The transient electrochemical response to a large dc 
voltage (without Faradaic reactions) has only been ana- 
lyzed quite recently [l], [52|, [53|, [54], [5f| ■ Even in the limit 
of thin double layers, the large voltage leads to a number 
of new dynamical effects not captured by circuit models. 
Additional time scales enter the problem, other than the 
fundamental RC time scale (which can be expressed as 
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XdL/D where D is the ion diffusivity In the sim- 
plest one-dimensional problem with parallel-plate, block- 
ing electrodes, over-charging of the double-layer "capac- 
itors" leads to net adsorption of neutral salt from the 
bulk,regardless of the polarity of the diffuse-layer volt- 
age PJ. This process is coupled to slow diffusive relax- 
ation of the bulk concentration (at the time scale L 2 /D), 
which leads to transient concentration polarization and 
thus breakdown of Ohm's law for the bulk "resistor" . 
In higher dimensions, large applied voltages also trigger 
surface transport of ions through the double layers [56| , 
which completes flux loops driven by bulk concentration 
gradients in and out of the double layers [52| . 

At large voltages, another important consideration is 
the breakdown of dilute-solution theory [H?], HI] , includ- 
ing the Poisson Boltzmann (PB) model of the double 
layer [53[, and, more generally, the PNP equations from 
which it is derived [H[ . These classical models are strictly 
valid only for a dilute solution of point-like ions, but, 
even in a very dilute bulk solution, the application of 
a large voltage can lead to the crowding of counterions 
near a highly charged electrode. Among many possible 
modifications of the PB model, one must at least ac- 
count for the finite sizes of ions and solvent molecules. 
This generally leads to the formation of a condensed layer 
of crowded ions, anticipated by Stern [59| and first de- 
scribed 1942 by Bikerman [6(|, whose simple modified 
PB (MPB) model has been rederived several times in 
different contexts [U|, H3, H [H, HE J6(| [67[ . As per- 
haps first predicted by Freise in 1952 (68|, the widening 
of the condensed layer generally causes the diffuse-layer 
differential capacitance to decay at large voltages - the 
opposite trend from PB theory, which allows ions to pile 
up with exponentially diverging concentration. This has 
major implications for the dynamics of electrolytes at 
large voltage s I53L 1571 [5a. 16911 as well as ionic liquids and 
molten salts [7fJ, lLU, Um (where crowding dominates 
in the absence of a solvent). In electrolytes, for the same 
reason, steric constraints also greatly reduce salt adsorp- 
tion and surface conduction compared to PB theor y by 
limiting the charge density of the double layer (53l . |56{ . 
All of these conclusions are independent of the model for 
steric effects on the chemical potential of ions in a con- 
centrated solution [57], [z3| and can be extended to more 
general situations, without assuming thin double layers, 
by deriving modified PNP (MPNP) equations [H[. 

A number of recent developments provide further mo- 
tivation for our work. In a recent paper [55|, Beunis 
et al. have revisited the problem of a suddenly applied 
large dc voltage in a blocking cell and studied the forma- 
tion of transient space-charge layers at very large volt- 
age, a possibility predicted in Ref. [l[ and analyzed pre- 
liminarily in Ref. [75] ■ Two recent papers, by Suh and 
Kang [jg, [77[ , analyze the weakly nonlinear response of 
an electrolyte to an ac voltage, which is relevant for 
many of the experimental situations described above. 
By coupling weakly nonlinear, charge relaxation to fluid 
flow, novel concentrated-solution effects can enter the- 



ory of induced-charge electrokinetic phenomena in large 
ac voltages [Hj]. In the context of electrodialysis mem- 
branes, it is well known that strongly nonlinear effects 
are important and can lead to electro-osmotic instability 
at the limiting current 0, [5(| [5lj . but this possibil- 
ity is just beginning to be explored experimentally us- 
ing large ac voltages. Building on recent observations of 
salt depletion and electro-convection near micro/nano- 
channel junctions [78| . the Rubinstein-Zaltzman insta- 
bility has been been demonstrated experimentally by ap- 
plying low-frequency ac (square-wave) voltages to con- 
fine it to slowly oscillating boundary layers j79[. This 
experiment raises interesting theoretical questions about 
the periodic breakdown and restoration of the quasi- 
cquilibrium structure of the double-layer under strong 
ac forcing, which are a major focus of this paper. 

In this work, we analyze the strongly nonlinear, time- 
dependent response of an electrolyte or ionic liquid to 
a large ac voltage, apparently for the first time. The 
imposition of a time scale (the ac period) is a significant 
complication compared to case of a sudden dc voltage, so 
we focus on the simplest geometry of parallel-plate block- 
ing electrodes and ignore any transverse instability. Fol- 
lowing Ref. [l| , we analyze the resulting one-dimensional 
problem starting from the classical PNP equations and 
derive accurate asymptotic approximations for thin dou- 
ble layers. We also consider the MPNP equations of 
Ref. J5J] to highlight steric effects under ac forcing. Using 
both PNP and MPNP models, we study the formation 
and collapse of transient space-charge layers at large volt- 
ages. While Beunis et.al [55| focus on the extreme case 
where the space-charge layer completely dominates the 
response at very large voltages (in sufficiently large sys- 
tems and high salt concentrations), we aim to derive a 
reduced model that is uniformly valid for all voltages and 
all salt concentrations, ranging from dilute electrolytes to 
concentrated solutions and ionic liquids. In spite of the 
mathematical complexity of these problems, our goal is 
to extract generic predictions and useful analytical ap- 
proximations to aid in interpreting experimental data. 

The paper is organized as follows. We begin in Sec. [U 
by stating the mathematical problem, converting to di- 
mensionless form, and showing full numerical solutions 
used to test our subsequent analytical approximations. 
In Sec. IIIII we briefly go though the asymptotic analy- 
sis for double layers in quasi-equilibrium, adapting the 
results of Ref. [lj concerning the transient dynamics, to 
our case of interest, namely the steady-state response 
when an ac voltage with frequency around the inverse 
RC time is applied. In Sees. [IV] and [V] we study the 
dynamic response in the weakly and strongly nonlinear 
regimes, respectively, and discuss how the circuit model 
is changed when we go from the weakly to the strongly 
nonlinear regimes. We also compare the strongly nonlin- 
ear asymptotic analysis to the full numerical solution. In 
Sec. I VII we develop an asymptotic analysis for the case 
when the double layers are driven out of quasi-equilibrium 
to form bulk space-charge and also compare those re- 
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suits to the full numerical solution. Finally, in Sec. IVIII 
we summarize and briefly discuss extensions to higher 
dimensions, Faradaic currents, and nonlinear electroos- 
motic flows, building on the initial study of Ref. [75|, 
and we leave the reader with some open questions. 



II. GOVERNING EQUATIONS 

A. General models 

In any continuum model, the transport of ions in the 
electrolyte is governed by a mass conservation law 



8 tCl = -V • Pi, 



(1) 



where q is the local concentration of the ith ionic species, 
Fj is the flux, and we neglect any bulk reactions in the 
electrolyte which could produce or consume ions. Quite 
generally, in a concentrated solution, the flux can be ex- 
pressed in terms of the gradient of the electrochemical 
potential \ii as 



A 

kT 



UCi 



(2) 



where the first term describes ion transport by diffu- 
sion and electromigration with diffusivity Di, k is Boltz- 
mann's constant, T is absolute temperature, and the sec- 
ond term describes advection at the mean fluid velocity 
u, as determined by momentum conservation. For a di- 
lute solution, the chemical potential fa takes the ideal 
form, with contributions from entropy and mean electro- 
static energy, 



Mi 



kT log ^ + Zieq 



(3) 



where <fi is the electrostatic potential, 2j is the ionic va- 
lence and e the electron charge. Equation ((T|) then re- 
duces to the Nernst-Planck equation 



-D 



(V, 



kT 



(4) 



In the usual mean-field approximation, the electrostatic 
potential is self-consistently determined by the charge 
density p through Poisson's equation 



V • (eV<; 



E 



2 / C Cj 



(5) 



where e is the electrolyte permittivity, which we take to 
be constant. This completes the classical PNP equations, 
which underly most of electrochemical transport theory. 
As noted above, the characteristic length scale in these 
equations (the Debye-Hiickel screening length) is 



ekT 



(6) 



where c* is the nominal bulk concentration of the ith 
ionic species. 

In the present work we focus on dilute electrolytes 
for which the nominal bulk salt concentration is small, 
seemingly within the range of applicability of the PNP 
equations. Even in a very dilute bulk solution, however, 
when a large external bias is placed on the electrodes in 
the system (only a few times kT/e), ions accumulate at 
the surface, and the dilute-solution approximation must 
break down [H, H3, HI] ■ Following Kilic et al. [HI], we 
will solve modified (MPNP) equations based on the old- 
est and simplest approach to steric effects of ion crowd- 
ing of Bikerman [60( , which corresponds to the following 
model for the chemical potential in a binary z : z elec- 
trolyte [Hlzl, 



/i ± = fcTlogc± ± zecj) — fcTlog(l — 



c_a 3 ), (7) 



where a is an effective molecular length scale. (For a 
history of this model and related concentrated-solution 
theories, see Ref. [58(.) The correction term, which can 
be interpreted as an activity coefficient fa = exp[(/Xj — 
/^ deal )/fcT], is related to the entropy of the solvent 
molecules and imposes a maximum ion concentration 
Cmax = a~ 3 ; it can be derived from the statistical me- 
chanics of of equal-sized ions and solvent molecules on a 
cubic lattice of spacing a in the continuum limit. In equi- 
librium, fi± = constant, the ions effectively obey Fermi- 
Dirac statistics, rather than classical Boltzmann statis- 
tics, due to the excluded volume effect [gj [62|, [63|, [H, 
[H, 0, [67L f70ll . Although more sophisticated models for 
fii exist 54 , 5^, [74[ , this approach at least qualitatively 



captures the effects of volume constraints with only one 
additional parameter a. 

For boundary conditions at the (blocking) electrodes, 
we assume no electrochemical reactions, so the normal 
ionic fluxes must vanish n ■ Fj = 0. To close the system, 
we follow many prior authors [E [M [Q, [ij, |H, M, HI 
and allow for a compact (Stern) layer or thin dielectric 
coating separating the electrode from the electrolyte with 
a constant "surface capacitance" per unit area Cg, which 
leads to a mixed boundary condition 



C s (V cxt -0)+en-V<j ) = O. 



(8) 



Here n is a surface normal pointing into the electrolyte, 
Cs = £s/hs can be ascribed to a surface coating of thick- 
ness h$ and dielectric constant eg, and V ex ±(t) is the ex- 
ternal potential applied at the electrode. 

For the present analysis we focus on a symmetric bi- 
nary electrolyte with equal diffusivity D + = D— = D 
and valence z + — — z_ = z for the two ionic species. 
Moreover we restrict our attention to the simplest pro- 
totypical microelectrochemical system, consisting of the 
electrolyte confined between two parallel, planar, block- 
ing electrodes at x = ±L, as sketched in Fig. [T] By 
symmetry, this rules out any effects of surface conduc- 
52l . l5r| or ac electroosmotic flow 0, [l(| E| and al- 
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lows us to focus on the strongly nonlinear response due to 
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FIG. 1: Sketch of ID model problem. The electrolyte 
is confined between parallel-plate blocking electrodes sepa- 
rated by a gap of width 2L, and a harmonic potential of 
V e xt(t) = ±V sin(uit) is applied to the left and right elec- 
trode, respectively, so the overall potential drop across the 
cell is 2V sin(wt); this corresponds to a 4V peak-to-peak volt- 
age or mis. 

the excessive accumulation of ions in the screening layers 
at the electrodes. In summary, the system is identical to 
that studied in Refs. Q (PNP) and [Hi] (MPNP), except 
that we apply an ac voltage rather than a step dc voltage 
and study the periodic response after all transients have 
decayed. We shall see that imposing an external time 
scale (the ac period) fundamentally alters the dynamics 
and complicates the analysis. 

B. Dimensionless form in one dimension 

We cast the problem into dimensionless form using L 
as the reference length scale and the RC relaxation time 
t = X]jL/D as reference time scale [jj, so that time and 
space are represented by the dimensionless variables t' = 
t/r and x' = x/L. The potential and ionic concentrations 
are rescaled as <fi' = 4>ze/kT and c' ± = c±/c*, where kT/e 
is the thermal voltage scale and c* is the nominal bulk 
electrolyte concentration. 

After dropping the primes from the dimensionless vari- 
ables, the governing equations take the form 

-e 2 c^=^(c+- C -), (9) 

d t c± = -ed x F ± , (10) 

where the fluxes F± are given by 

F± = -c±d x p±. (11) 

For a dilute electrolyte the electrochemical potentials re- 
duce to 

p± = log c±±<p, (12) 

and we arrive at the Nernst-Planck equations in dimen- 
sionless form 

d t c± = ed x (d x c± ± c ± d x (f>). (13) 



When steric exclusion is taken into account we get 

p± = log c± ± (j) - log(l - uc), (14) 

where the parameter u = 2c* a 3 is the nominal volume 
fraction of the ions in the electrolyte [53| . 

It is convenient to introduce also the average ion or 
"salt" concentration and (half) the charge density 

c=i(c+ + c_) and p=-(c+-c_), (15) 

in terms of which the transport equations can be rewrit- 
ten as 

d t c = -ed x F, (16) 
d tP = —ed x J. (17) 

Here F = \{F + + F_) and J = \{F+ - FJ) are the 
average salt flux and current density, respectively, 

F = -d x c/(l - uc) - pd x ct>, (18) 
J = -dxp - cd x 4> - vpd x c/(l - vc). (19) 

Since we assume blocking electrodes we have no-flux 
boundary condition at the electrodes, 

F± = or F = J = 0, (20) 

whereas the compact layer b.c. reduces to 

V cxt -4>= T^5d x <j) at x = ±1. (21) 

Here V ey ±(t) = ^fVsmicvt) is the electrode potential and 
6 = C'd/Cs is the ratio of the compact-layer capacitance 
Cs = £s/hs to that of the diffuse layer Cd = £/A_d in 
the low voltage limit. 

C. Dimensionless parameters 

The PNP model contains three dimensionless parame- 
ters: e, V, and S. In aqueous electrolytes, the screening 
length has submicron scale, Ad ~ 1 — 100 nm, so the 
diffuse-charge boundary layers near the electrodes typi- 
cally have a very small dimensionless width e = / L <C 
1, at least in microsystems where L > 1 pm.. Although 
e > 1 is possible in nanosystems, we restrict our atten- 
tion to the typical case e 1, which is the basis for our 
asymptotic analysis. 

Contrary to most prior work, we focus on the nonlinear 
regime of large applied voltages, V » 1 (or, with units, 
V > kT/e w 25 mV), as in Refs. Q, IH, gl H3, [H, H . 
Since applied voltages larger than a few volts tend to 
trigger Faradaic reactions in aqueous electrolytes at ac 
frequencies around the inverse RC time [l2l [83 | , we en- 
vision experimentally relevant values of V f=a 1 — 200, 
although larger voltages can be sustained at higher fre- 
quencies or in non-aqueous solvents or liquid salts. Unlike 
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prior work, we allow for large enough voltages that the 
double-layers lose their quasi-equilibrium structure. 

The parameter 5 = ehs/es^D = ^s/^d (where A5 = 
hse/es is an effective thickness for the compact layer) 
can be estimated in some cases, but it is usually adjusted 
to fit experimental data. For example, let us consider 
different surfaces in contact with a 1 mM aqueous elec- 
trolyte with Ad = 10 nm. For a thin dielectic coating, 
such as a natural TiC>2 oxide layer with hs = 4 nm and 
es = 110 where a constant Cs seems reasonable, we get 
S w 0.3, although much thicker dielectric la yer s yielding 

5 « 10 can arise in patterned microsystems [17] . In non- 
linear electrokinetic phenomena, the inferred value of S, 
required to match the standard dilute-solution model to 
experimental data, can be up to several orders of mag- 
nitude larger, although this is more likely due to failures 
of the model, and not directly related to surface capaci- 
tance [581 ] . 

With a dielectric electrode coating, the surface capac- 
itance is relativel y c lear, but in the classical picture of 
the Stern model [59|], it is associated with a hypothet- 
ical flat monolayer of water molecules, which limit the 
approach of hydrated ions at the metal/electrolyte inter- 
face. In that case, one would expect hs 1 A, and, since 
alignment of water dipoles is assumed to reduce the per- 
mittivity in the Stern layer to £5 as O.le 83, [Hj], we esti- 
mate S fts 0.1. The Stern picture, however, is complicated 
(at least) by electronic boundary layers in the metal [85j ] . 
chemisorption from the solution [86|, nanoscale surface 
roughness [13, HH, and crowding effects absent in the 
PNP model of the diffuse layer, all of which can be mis- 
attributed to the Stern layer, as emphasized by Bazant 
et al. [58]. Even for a smooth liquid-mercury electrode, 
the inferred Stern layer capacitance is voltage depen- 
dent [H, (86|. Nevertheless, since our goal here is to 
analyze the nonlinear dynamics of ions in solution, we 
will simply assume a constant surface capacitance and 
allow for a wide range of values 6 = 0.01 — 10. 

In the simple MPNP model there is one more dimen- 
sionless parameter, v = 2a 3 c* , which controls the im- 
portance of crowding effects. The nominal concentration 
c* could range from a very dilute 1 /iM solution with 

6 x 10 20 ions/m 3 to a concentrated 1 M solution (near 
physiological salt levels) with 6 x 10 26 ions/m 3 . A natu- 
ral choice for the effective molecular lattice spacing a is 
the diameter of a hydrated ion, around 4 — 5 A for small 
ions in water, which would yield v w 10~ 7 — 10 _1 . Tak- 
ing into account the under-estimation of steric effects in 
a hard-sphere liquid by our lattice-based model [H, [74[ , 
the value of a could be increased by roughly a factor of 
two [6{|. Electrostatic correlations also become impor- 
tant when ions are crowded at this scale, comparable to 
the Bjerrum length, 7 A in bulk water. As a crude ap- 
proximation, therefore, we may consider v as large as 0.4 
in a concentrated electrolyte. 




FIG. 2: Numerical solution of the PNP eqs. for V = 30, 
u = 0.3, 8 — 0.3, and e = 0.001. (a) Potential variation in 
time and space; solid black line shows the external potential 
on the electrodes V cyL t. The inset zooms onto the rapid poten- 
tial variation across the diffuse screening layer, (b) Zoom on 
cation concentration near the electrodes; anion concentration 
is identical but phase shifted one half period in time. 



D. Numerical solution 

Before embarking on our asymptotic analysis, we 
present some numerical solutions of the PNP model for 
the problem sketched in Fig.[TJ which will be used to test 
and calibrate various analytical approximations below. 
As noted above, steric effects in the MPNP model tend 
to reduce nonlinearities, so the PNP model serves as a 
more stringent test case. 

We use the Comsol finite element package [8!| to solve 
numerically the PNP model in the form of Eqs. (J5J), l|16[). 
and {TTJ) with boundary conditions (J20J) and (J2TJ). It 
is necessary to use a very fine mesh of Ax m i n ~ 10~ J 
close to the electrodes in order to resolve the highly com- 
pressed (and unphysical) diffuse-layer structure in the 
PNP model, even at e = 10~ 3 . For our ID problem, this 
is straightforward to achieve using a nonuniform graded 
mesh, but in 2D or 3D it would pose a serious problem. In 
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FIG. 3: Numerical solution of the PNP eqs. for V = 120, 
(j = 0.3, 8 — 0.3, and e = 0.001. (a) Potential variation in 
time and space, (b) Zoom on cation concentration near the 
electrodes. 



the ion concentration varies very rapidly from a max- 
imum of max{c±} « 3 x 10 3 at the electrode surface 
down to around unity over a distance of O(e). The inset 
in Fig. Ufa) zooms onto the rapid potential variation in 
the screening layer, and the visible difference between the 
potential on the electrodes and in the electrolyte corre- 
sponds to the compact layer voltage, cf. Eq. |2ip. 

Within a distance of about 0.1 from the electrodes we 
see a nonuniform pattern in the cation concentration pro- 
file that oscillates at twice the driving frequency, 2u>. 
The concentration has a minimum just about the time 
when the screening layer is fully charged, and a maxi- 
mum when the screening layer changes polarity (occurs 
for t sa 1.75 and again at 12.25). The anion concentra- 
tion shows a fully similar pattern so that effectively this 
"diffusion layer" is charge neutral. 

Figure[3]shows the solution for V = 120 with otherwise 
the same parameters as in Fig. [5J The overall picture is 
essentially the same as before: The bulk ionic concentra- 
tions are constant, but due to the massive accumulation 
of ions around the electrodes (maximal concentration ex- 
ceeds 4 x 10 4 ), the bulk concentration is down to 0.86, 
i.e., 14% below the nominal value. 

An interesting feature is seen in Fig.[3fb) for t w 5.24: 
Close to the electrode at x = — 1 there is an extended re- 
gion where the cation concentration drops to zero, while 
at the same time the anion concentration is also low but 
clearly nonzero (see Fig. [TJ] below for a more detailed 
view). This transient "space-charge layer" is similar 
to the steady counterpart described by Rubinstein and 
Shtilman for the case of dc Faradaic conduction when an 
electrochemical cell is driven above the diffusion limited 
current [47],[48| . It is also clear from Fig.[3]Ja) that there is 
a significant potential drop across the space-charge layer. 



fact, overcoming such limitations is a major motivation 
for our development of accurate boundary-layer approxi- 
mations below. The steady-state periodic response is ob- 
tained by integrating forward in time, using the default 
time-dependent solver of Comsol. Since the transient 
diffusive relaxation in the bulk is slow, it is necessary to 
integrate for a very long time, up to 100 times the period 
of the driving voltage or more, before the steady-state 
periodic solution is reached. 

Figure [5] shows the result for V = 30, co = 0.3, S = 0.3, 
and e = 0.001: Fig. [2fa) shows the potential <j>(x,t), 
and Fig. Hfb) displays the cation concentration profile 
C-\-(x,t). In the bulk region both the cation and anion 
concentrations are constant and (very close to) unity, 
and the electrolyte therefore behaves like an ideal resis- 
tive medium with unit conductivity. The potential varies 
linearly throughout the bulk region, driving a constant 
ohmic current density, and shows a roughly harmonic 
time variation that is about 45° ahead of the external 
potential V cxt = ±Vsin(u;i) (solid black line at x = =pl). 

In the diffuse screening layer close to the electrodes 



III. ASYMPTOTIC ANALYSIS 

A. Nested boundary layers 

In the limit of thin double layers, e = Afl/L < 1, the 
dynamical problem can be analysed by matched asymp- 
totic expansions QJ. The standard procedure begins by 
seeking regular expansions in the form of power series 

c = c (0) +ec (D +e 2 c ( 2 )+..., (22) 

substituting into the governing equations, and collecting 
like powers of e. This procedure is guaranteed to converge 
in the limit e — > with all other parameters held fixed. 
However, for any fixed e > there could be e-dependent 
restrictions on the other parameters, in particular the 
driving voltage V , for a truncated expansion to produce 
accurate results. Following Bazant et al. [lj, 52] we denote 
the regime where such conditions hold as "weakly non- 
linear" , as opposed to the "strongly nonlinear" regime 
where the standard asymptotic expansions breaks down. 
It is on this "strongly nonlinear" regime that we focus 
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our attention. We aim at deriving the leading-order dom- 
inant balance in the joint limit e -> and V — ► oo, and 
since we focus exclusively on the leading-order approx- 
imation, we drop the superscript ^ on all variables in 
order to simplify the notation. 

For small applied voltages, it is well known that in the 
limit e — ► the diffuse part of the double layer acts as 
a mathematical boundary layer of O(l) non-zero charge 
density and 0{e) thickness on the leading-order quasi- 
electroneutral bulk region at the O(l) length scale of the 
geometry. This is the mathematical justification for lin- 
ear circuit models. In the case of blocking electrodes, the 
characteristic RC time scale for charging of the double 
layers is O(l) in our dimcnsionless units. 

The application of a large voltage leads to the new ef- 
fect of salt adsorption by the diffuse layer and related 
depletion of the bulk concentration, first described by 
Bazant et al. and in higher dimensions, surface con- 
duction through the diffuse layer becomes important at 
the same time [52| . For a suddenly applied dc voltage at 
blocking electrodes, during the initial RC charging phase 
over O(l) time, a thin quasi-electroneutral diffusion layer 
extends to 0{\ft) width. The next phase of relaxation 
proceeds at the slow 0(e _1 ) time scale for bulk diffusion, 
as concentration gradients spread across the cell to 0(1) 
distances. 

For a large applied ac voltage, this picture is altered 
by the imposed time scale. In the case of ac forcing close 
the RC time scale, uj^ 1 = 0(1), the oscillating voltage 
generally leads to the formation of a thin, nested oscil- 
lating diffusion layer confined to steady 0(^/e) thickness. 
Since salt adsorption by the inner diffuse layer is positive, 
regardless of the polarity, the diffusion layer oscillates at 
twice the driving frequency. It is also accompanied by 
a gradual bulk salt depletion that propagates across the 
cell over O(l) times after the ac voltage is turned on, 
similar to the case of the sudden dc voltage. However, 
in this work, we ignore such initial transients and focus 
on the steady ac response after the bulk has relaxed to 
steady state. 

In this way, we are lead to analyze a new, nested 
boundary-layer structure sketched in Figure 21 consisting 
of the "outer" bulk region (unit length scale) , a "middle" 
diffusion layer (i/e length scale), and the "inner" diffuse 
part of the double layer (e length scale). This picture 
remains valid until the voltage becomes large enough to 
fully deplete the middle diffusion layer, leading to the 
formation of transient space-charge layers extending by 
0(y/e) or more into the cell, twice per ac period. Our 
goal in the rest of the paper is to develop uniformly valid 
asymptotic boundary-layer approximations in all of these 
cases. For clarity, we denote asymptotic approximations 
for each region by different accents, as indicated Figure 
2J For example, the salt concentration c is asymptotic 
to c in the bulk, c in the diffusion layer, c in the diffuse 
charge layer, and c in the space-charge layer. 





- "Outer" bulk electrolyte 
c, J 




_ "Middle" diffusion layer 


Space-charge layer 
c, $, Vo 




. "Inner" diffuse layer 

^' ^' e tc. f Condensed layer u 


Electrode Compact layer S 



FIG. 4: (Color online) Schematic picture of nested boundary 
layers in matched asymptotic expansion: The bulk "outer" 
region is connected via the "middle" diffusion layer to the 
"inner" diffuse layer. A compact (Stern) layer separates the 
electrolyte from the blocking electrode. At large voltage, lo- 
cal salt depletion in the diffusion layer can cause the double 
layer to change to a non-equilibrium structure, with an ex- 
tended "space-charge" layer that is completely depleted of 
coions. Further, when the concentration in the diffuse layer 
approaches the steric limit, a condensed phase of ions forms 
at the electrode. The figure also indicates some of the vari- 
ables and parameters introduced in the asymptotic analysis 
of the different layers. 



B. Quasi-electroneutral bulk 

We begin by analyzing the solution in the bulk region. 
The Poisson equation ([9]) shows that the charge density 
vanishes to both zeroth and first order in e so that at 
leading order 

c + = c_ = c, (23) 

denoting bulk variables by a bar accent. The bulk salt 
concentration displays diffusive dynamics on the time 
scale t — et [l[, but on the RC time scale the concen- 
tration profile is constant in time c = c(x). Here we 
focus on the steady state after the bulk transients have 
relaxed, and by symmetry of our simple model problem, 
the bulk concentration is then simply constant in space, 
c = c Q . The leading order potential varies linearly in 
space 

$ = -^-x, (24) 

Co 

where J(t) is the ohmic current and c acts as the bulk 
conductivity. 

In the weakly nonlinear regime the bulk concentration 
is at the nominal value, c Q — 1, whereas in the strongly 
nonlinear regime the adsorption of ions in the double 
layers may be so strong as to induce c < 1. 

The problem has the following symmetries about the 
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origin 



4>(x,t) 
p(x,t) 
c(x, t) 




(25) 



In Sections IIII CIIIII El below we focus on the nested 
boundary layers developing at the left electrode, and for 
convenience we therefore perform a change of variables, 
y = 1 + x, such that y — corresponds to the electrode 
surface and y > to the interior of the cell. 



and we obtain the standard PB equation 

<9~-0 = c s sinh?/>, (31) 
yielding the familiar Gouy-Chapman (GC) solution 

-l 



ip = 4 tanh 



tanh(C/4) 



,-vc. y 



(32) 



Here the integration constant ( = i(j(0) is simply the lead- 
ing order zeta potential, and 1 / \ft s is the local effective 
Debye length. 



C. Quasi-equilibrium double layer 

The singular perturbation in the Poisson equation (|9|) 
gives rise to a boundary layer of width O(e) where the 
charge density is nonzero to zeroth order in e. Intro- 
ducing a scaled spatial variable y = y/e to remove the 
singular perturbation, we can seek regular asymptotic 
expansions (denoted by tilde accents) in the "inner" dif- 
fuse layer. Substituting into Eqs. (fT0|) and (fTTj) and us- 
ing dy = ed y , we find that the double layer is in quasi- 
equilibrium at leading order with constant electrochem- 
ical potential fj,± across it. The value of jx± is deter- 
mined by matching with the solution in the adjacent 
quasi-electroneutral diffusion layer 



lim {log c ± (f)}. 



(26) 



The quasi-equilibrium arises because the diffuse charge 
dynamics relaxes on the Debye time scale t = t/e, which 
is much faster than the RC charging time. The ion dis- 
tributions are determined from Eq. (| 14[) as 



1 + i/c s (coship — 1) 



(27) 



where ip — (f> — cj> is the excess potential in the double 
layer relative to the diffusion layer, and c s is the limiting 
value of the salt concentration c s — lim^o c as seen from 
the double layer. The excess potential satisfies the MPB 
equation 



c s sinh ip 



1 + uc s (coshtp — 1) 
which can be integrated once to get the field [H, \& 



(28) 



dy^P = -sign(V0y 21og[l + 2t/c s sinh^(7/;/2)]/^. (29) 

Note that the ion concentrations in Eq. (|27]) are bounded 
above by steric exclusion, c± < 2/u, while at much lower 
concentrations they reduce to the usual results from di- 
lute theory: In this limit (y — -> 0) the ion profiles are 
given by the Boltzmann equilibrium distribution, 



c s e 



TV 



(30) 



D. Surface conservation laws 

The redistribution of ions across the diffuse layer is 
instantaneous on the RC time scale, but the total amount 
of ions absorbed can change only by flux into the layer 
from the adjacent diffusion layer. Following Bazant et 
al. [H, [53, [13, [El] we quantify this by considering the 
excess amount of each ionic species accumulated in the 
double layer, w± — ew±, where 



1 



w± = 



/ (c± - c±) dy = (c± - c±) dy. (33) 

Jd.l. Jo 



The time evolution of w± is then determined by 



dtw± = I d t (c± - c±) dy = - lim F± (34) 
/o v^°° 



lim F± , 



(35) 



where the last equality is obtained by flux matching be- 
tween the double layer and diffusion layer. We also define 
the diffuse charge and excess salt concentration by 



1 



(w + — W-) and 



1 



(w++w^). (36) 



Since the diffusion layer is quasi-electroneutral at lead- 
ing order, the double-layer charging process is coupled 
directly to the bulk electric current 



d t q = -J(t). 



(37) 



Variations in the excess salt w are coupled to the dynam- 
ics in the diffusion layer 



d t w 



lim F 



-F n 



(38) 



where the flux injection F a {t) at the inner "edge" of the 
diffusion layer should be understood as the driving force 
behind the oscillations in the salt concentration in the dif- 
fusion layer. These relations exemplify the general math- 
ematical theory of surface conservation laws, in which 
the total excess concentrations in a diffuse interface are 
coupled to normal (and surface) fluxes in a concentrated 
solution 15611. 
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E. Oscillating diffusion layer 

The diffuse screening layers at the electrodes period- 
ically absorb and expel an excess amount of ions from 
the surrounding electrolyte. However, the bulk trans- 
port of neutral salt is essentially a diffusion process on 
the time scale t — ei, which is much slower than the ac 
driving that we consider here, and hence the leading or- 
der dynamics are confined to a diffusion layer of 0(- l /e) 
width around the electrode [l|, • We therefore intro- 
duce a scaled spatial variable y = y/yfe and seek regular 
asymptotic expansions (denoted by hat accents) in this 
"middle" diffusion layer. Substituting into the Poisson 
equation ([9]) we find that, like in the bulk, the charge 
density vanishes to zeroth order in e, so that the leading 
order ion concentrations are equal 



c+ 



(39) 



Equation (| 1 6|) then reduces to a simple diffusion problem 

d t c = dp, (40) 

to be solved on the interval y £ [0,oo). Matching to 
the bulk solution requires lim^oo c ~ lim^o c = c D , 
whereas the boundary condition at the inner edge of the 
diffusion layer is determined by matching with the salt 
flux out of the double layer, cf. Eq. 



(41) 



—= lim d{,c = lim F = F . 



The solution to the ID diffusion problem can be ex- 
pressed in terms of a convolution integral [l[ 



where 



G(y,t) 



G(y,t-t')F (t')dt', 



-y 2 /4t 



nt 



(42) 



(43) 



is the Green's function for the diffusion equation with a 
sudden unit flux at t = + injected at the boundary, 



G(y,0) = 0, -%G(0+,i) = S + (t). 



(44) 



for a semi-infinite domain. Technically, Eq. (|4"3")l is the 
first term in an expansion for the Green function in a 
finite bulk domain (Eq. 24 of Ref. [1]), which would 
be needed to describe the initial transient when the ac 
voltage is first turned on. Here, we focus on the steady- 
state response, after initial diffusion layers have relaxed 
across the cell, thereby lowering the uniform bulk concen- 
tration c (see below), and the oscillating diffusion layers 
have only 0(\fe) width, consistent with the semi-infinite 
approximation (|43[) . 

To describe this situation, since the flux injection is 
periodic we may rewrite Eq. (|42[) as 



(45) 



c = c + v^/ G u {y,t-t')F {t')dt', 



where T = 2-k/lo is the driving period and 
1 



T 



(46) 



n=l ^ 7 

is the Green's function for a periodic influx of salt, 

oo 

(G u (0,i)) = 0, -a t G u (Q+,t)= S+(t-nT). (47) 



n— — oo 



Equations Eqs. (|4"2"1) or (I45|) clearly show that in the 
weakly nonlinear regime, where F Q is 0(1), the concen- 
tration in the diffusion layer c is equal to the bulk c at 
leading order; the flux injection only gives rise to an 
0(-y/e) perturbation. The strongly nonlinear regime is 
essentially defined as the regime of driving voltages high 
enough that w and F Q grows to 0(1/ y/e) and the varia- 
tions in c reach O(l). 

Since the diffusion layer is charge neutral at leading 
order, the current is constant across it and equal to the 
bulk current J(t). However, the conductivity differs from 
its bulk value, which gives rise to transient concentration 
polarization. There is an excess electrostatic potential 
variation ip = <f> — (f>, and an excess field given by 



1 - -fl 
-jfy^ = J[ 



(48) 



In the weakly nonlinear regime when (1/c— 1/c) is O(^fe) 
it is clear that tp is only an 0(e) perturbation; in the 
strongly nonlinear regime ip grows to 0(y/eJ) perturba- 
tion which is, however, still negligible compared to the 
bulk 4> = —Jx/c . 

Finally, the leading order charge density in the diffu- 
sion layer can be evaluated by substituting Eq. (I48[) into 
the Poisson equation to get 



P : 



3/2 JdyC 



(49) 



The quasi-electroneutral solution in the diffusion layer 
remains valid for \p\ <C c; we return to this aspect in 

sec. rvn 



F. Closing the problem 

In order to close the coupled problem for the dynamical 
variables J(t), q(t), ((t), w(t), and c s (t) we need a few 
more relations between them. The first is obtained by 
writing the overall potential drop over the boundary lay- 
ers, from the electrode to the bulk electrolyte at x = — 1, 
as the sum of the contributions from the compact and 
diffuse layers, to get 



-l,t) = V cxt - J/c = -q6 + (. 



(50) 



Since we focus on the leading order approximation, we 
neglect here the small potential drop over the diffusion 
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layer. Next, the diffuse-layer voltage C can be related to 
the diffuse charge through Eq. (|2T)|) for the field at the 
electrode surface, yielding 

q = -sign(C) ^21og[l + 2i/c a sinh 2 (C/2)]/z/, (51) 
and in the dilute limit this reduces to Chapman's formula 
q = -2 v / ^sinh(C/2). (52) 
The charge-voltage relation can be inverted to get 



C = -sign(g)2sinh- 1 \ • (53) 

The excess salt concentration can be expressed in integral 
form [HI] 

c - 

w= / -jd^ (54) 

and using PB theory in the dilute limit the integral can 
be evaluated to get [l| 

w = 4Vc^sinh 2 (C/4), (55) 

or, eliminating £ we obtain 

w = y/q 2 + 4c s - \/4ST- (56) 

For the MPB model, Eq. ([54]) is difficult to handle analyt- 
ically, but numerical integration shows that Eq. (I56p ap- 
proximates the integral well, with relative error of 0(v). 

The bulk salt concentration c is determined by im- 
posing the global conservation of salt in the cell 

J c(x,t)dx = 2. (57) 

As noted in Ref. [45j, integral constraints on the total 
number of inactive ions are generally required for steady- 
state problems to replace information about the initial 
condition (e.g. when the voltage is first turned on) that 
is preserved during time evolution with no-flux boundary 
conditions. 

There is a periodic exchange of salt between the inner 
diffuse and middle diffusion layers, but Eq. (|4"B")) shows 
that (c) = c, i.e., the diffusion layer does not contain any 
excess salt on time average. As a result, the uniform bulk 
concentration in the (time-periodic) steady state, c Q , is 
reduced only by the time-averaged salt adsorption of the 
diffuse layers, 

c = l-e{w). (58) 

which also describes the (static) steady state after a 
sudden dc voltage is imposed [l[. This, together with 
Eqs. ([37]), d3Hl) , (E01), and flU constitute a set 
of "ordinary" integro-differential-algebraic equations in 



time for the dynamical variables J(t), q(t), F D (t), w(t), 
c s (t), and C(i)- 

Our focus in the present work is on the steady-state 
periodic response, and we explicitly made use of this in 
deriving our dynamical model by replacing the transient 
Eq. ((431) with Eq. ([46]). The problem could be solved nu- 
merically by a relaxation method, representing each dy- 
namical variable by a truncated Fourier series, as done in 
Ref. [lH . Here, however, our approach is to integrate the 
dynamical equations by a timestepping algorithm; the 
integration is continued until a periodic state is reached, 
typically within 10-20 periods of the driving voltage. The 
timestepping approach is well suited for integrating also 
the non-equilibrium model developed in Sec. IVI1 and 
is much more efficient on computer memory for solving 
problems with 2D or 3D electrode geometry. Further de- 
tails are given in our supplementary material [90j . 



IV. WEAKLY NONLINEAR REGIME 

The "weakly nonlinear" regime defined in Ref. [l[ is 
characterized by fluctuations in the diffusion layer salt 
concentration being only a small perturbation to the bulk 
value, so that c = c = 1 at leading order. This is the re- 
sponse predicted by matched asymptotic expansions in 
the singular limit e — > with all other parameters held 
fixed, including V. As such the dimensionless, leading- 
order response is independent of e. We begin with an 
analysis of this regime, and the findings here form the 
basis for understanding the peculiarities of the "strongly 
nonlinear" regime in subsequent sections, where the so- 
lution has a nontrivial dependence on V and e. 



A. Charge-voltage relation 

The only nonlinearity in the weakly nonlinear model 
arises from the diffuse-layer charge-voltage relation, 
Eq. ([5Tj) . In Fig. [5] we plot the accumulated charge q 
as a function of the overall potential drop = ( — q5 
across the double layer for different values of the capaci- 
tance ratio 5 and nominal ion volume fraction v. 

In the Debye-Hiickel limit, £ <C 1, Eq. (p)Tj) can be lin- 
earized to get simply q = — C = — + 6). At larger 
voltage the classical PB theory predicts a dramatic in- 
crease in the diffuse-layer capacitance, and q grows expo- 
nentially with (. In Fig. [5] this behaviour shows directly 
on the curve 5 = v = (Gouy-Chapman (GC) model) 
that bends up sharply for ^ > 10; at "J = 20 the con- 
centration in the diffuse layer, cf. Eq. ([30]) . exceeds 10 8 
times the bulk concentration, which is absurdly high for 
aqueous electrolytes. 

This well-known unphysical artifact of PB theory is 
alleviated (but not eliminated) in the Gouy-Chapman- 
Stern (GCS) model by assuming a finite compact layer 
capacitance, corresponding to a positive value of S. Then 
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FIG. 5: (Color online) Quasi-steady (dimensionless) accumu- 
lated charge q in the double layer as a function of its voltage 
drop ^ = ( — q~8, plotted for different values of the Stern pa- 
rameter 8 with v = (solid) , and different values of the steric 
parameter v with 5 = (dashed), all in the weakly nonlinear 
regime with 6=1. 



at large voltage the major part is carried by the compact 
layer, * ss —qS, while the diffuse- layer voltage remains 
small, |C| ~ 2 log | g| « 21og I'F/JI. This regularizes the 
problem at moderate voltages, but the success may be 
misleading: It is unlikely that an angstrom-thick molec- 
ular Stern layer could withstand several volts without 
dielectric breakdown. Moreover, the GCS model does 
not impose a maximum charge density and at sufficiently 
large voltages still reaches unphysical ion concentrations. 

A more realistic approach could account for crowd- 
ing effects at large voltages using MPB theory 
e.g. leading to Bikerman's model described above 6C|. 



This is equivalent to PB theory at concentrations well 
below the limit of steric exclusion, c± <C 2/^, as seen 
clearly in Fig. [5l However, once steric effects saturate 
the charge density, the diffuse-layer capacitance quickly 
drops due to the condensed phase of ions forming at the 
electrode [H, [H [H, [H . This occurs for q 2 > 2/u, and 
at still larger voltage the overall potential drop is pri- 
marily on the condensed layer, with Eq. (fBTj) reducing 
to 




2 I*I * iti 2 
— — - for |*| > - 

v v 



(59) 



Comparing the GCS model for a Stern monolayer with 
S = 0.03 to Bikerman's model, the latter predicts (much) 
lower charging already for >F > 30 (or 750 mV) at v — 
10~ 4 . Even for our example of an oxide layer on the 
electrodes with 6 = 0.3, crowding effects in the liquid 
could significantly affect the charge-voltage response for 
"J > 100 (or 2.5 V), which is still within the range of 
many experiments. 

This asymptotic square-root dependence of the charge- 



voltage relation (|59)) is a generic consequence of volume 
constraints 58], not only in Bikerman's lattice-gas model, 
but also hard-sphere liquid models, since it corresponds 
to a diffuse layer of uniform charge density. Once the 
condensed layer forms, its voltage £ w vq 2 /2 can easily 
exceed that of the outer (PB) part of the diffuse layer 
Q sa 21og(q ; ), even while the latter remains thicker. In 
this regime of the model, the thickness of the condensed 
layer is t — £/e ~ vq, which does not become larger than 
the diffuse-layer thickness until q > \ jv. 



B. Dynamical response 

The leading order dynamic response in the weakly non- 
linear regime is governed by 



d t q = - J, 
-J = C-qS, 

C = 2sinlr 1 
c — c = 1. 



; ^ 2 /2 - 1 



2;/ 



(60) 
(61) 

(62) 
(63) 



This may be rewritten as a single ordinary differential 
equation for the double-layer voltage * [l[ 



Cd t <$ = J = V a 



(64) 



where C(^>) = — dq/dfy is the total differential capaci- 
tance of the double layer, and V ex t(t) = Vs'm(ujt) is the 
external driving voltage. 

We focus on the periodic response obtained by starting 
from an initially uncharged state and integrating forward 
in time until all transients have died out. Figure [H] shows 
the results for different values of the model parameters: 

Fig. [6{a) shows the result for V = 1, u> = 1, 5 = 
0.3, and v — 0. At this low voltage the charge-voltage 
relation is still essentially linear, so the system behaves 
like a linear RC circuit with time constant (l + S)^ 1 . The 
double-layer voltage * is dominated by the diffuse layer 
with the compact layer contributing only a small fraction 
8. 

Figure[Hb) shows the solution at larger voltage V = 30 
with u> = 0.5, 5 = 0.3, and v = 0. At this voltage the 
relation between £ and q is clearly nonlinear, £ stalls 
f° r |(?| si 10, and the double-layer voltage becomes dom- 
inated by the compact layer. When the double layer 
changes polarity this in turns makes the change of sign 
of ( look like a "sharp" transition which gives rise to a 
jump in the bulk current. Those features are even more 
pronounced in Fig. [6^c), showing the corresponding so- 
lution for (5 = 0, i.e., without any compact layer on the 
electrodes. The double-layer voltage remains low so the 
bulk current is almost in phase with the driving voltage. 

As discussed in the previous section, the very large 
capacitance of the diffuse layer predicted by PB theory 
is not realistic. For the solution in Fig. [SJc) the maxi- 
mal ion concentration in the diffuse layer almost reaches 
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FIG. 6: (Color online) Distribution of dimensionless cell voltage, 14 x t = J + C — Q$ (dotted), divided into contributions across 
the bulk electrolyte, J (solid), diffuse layer, £ (dashed), and compact layer, — qS (dash-dot), all scaled to the thermal voltage 
kT/ze, as a function of time for different values of the parameters V, u), 5, and v. The panels show (a) Debye-Hiickel limit, 
(b) Gouy-Chapman-Stern (GCS) model, (c) Gouy-Chapman (GC) model, and (d) Bikerman model. 



10 4 times the bulk concentration, which could easily trig- 
ger steric effects, even for a nominally dilute electrolyte. 
Figure [6^d) shows the result when such are taken into 
account with a bulk volume fraction v = 0.01. The re- 
sult is markedly different: When steric effects set in, the 
diffuse-layer capacitance drops and £ grows rapidly with 
q. At lower charging, though, the system is still governed 
by dilute theory, so we still see a rapid shift in £ with an 
associated jump in J when the double layer changes po- 
larity. 



C. Equivalent circuit 



A useful concept for analyzing the cell response is an 
equivalent circuit diagram like that shown in Fig. [7fa). 
The transport through the bulk electrolyte is represented 
by an ohmic resistor 2R = 2, and the charge accumula- 
tion in the double layer by a series coupling of two ca- 



(65) 



pacitors Cs and Cd 

III 

Here Cs — 1/6 is the capacitance of the compact (Stern) 
layer, and Cd = —dq/d(^ is the differential capacitance 
of the diffuse (Debye) layer, given by [H, [H, [68|, 170] ] 

| sinhCI 



C D = 

[1 + 2^sinh 2 (C/2)] ^2 log [l + 2v sinh 2 (C/2)] jv 

(66) 

In the Debye-Hiickel limit this reduces to Cd = I and 
C = 1/(1 + S). At higher voltage, PB theory predicts 
a dramatic increase of Cd — cosh(£/2), to the extent 
that C w 1/5. According to MPB theory, the diffuse- 
layer capacitance becomes a non-monotonic function of 
C, where the initial increase is followed by a decrease as 



Cd ~ l/vq_ ~ l/y2^C once steric exclusion sets in. 

The equivalent circuit representation is useful for un- 
derstanding and interpreting the system response. How- 
ever, from an experimental point of view the overall cell 
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FIG. 7: (Color online) (a) Equivalent circuit representation 
for weakly nonlinear dynamics: Compact and diffuse-layer 
capacitors in series with a bulk resistance, (b) Bode plot of the 
magnitude \Z\ and phase angle ZZ of the half-cell impedance 
for increasing driving voltage at S = 0.3 and v = 0. The 
characteristic frequency uj , where ZZ passes through —45° 
and \Z\ bends up, is marked with circles. 



impedance is a key property that can easily be measured 
with high accuracy, e.g., using a lock-in amplifier. We 
define the (half) cell impedance Z as the ratio between 
the first Fourier components of the applied voltage and 
the resulting current 



dt 



dt 



(67) 



Since the system is nonlinear, the impedance so defined 
is a function of both driving frequency and voltage. Fig- 
ure[7Jb) shows a Bode plot of the cell impedance Z for dif- 
ferent values of V at S — 0.3 and v = 0. The curve shape 
is characteristic of an RC series coupling. At high fre- 
quency the ohmic resistance of the bulk electrolyte dom- 
inates and \Z\ levels off at unity, while at low frequency 
the double-layer capacitance dominates and \Z\ oc uj^ 1 . 
At the same time the phase angle ZZ drops from zero 
at high frequency to —90° at low frequency. We define 
the characteristic frequency uj for a given driving voltage 
as that frequency where the phase angle passes through 
-45°, i.e., 




ZZ{u ) 



-45°. 



(68) 



FIG. 8: (Color online) Characteristic frequency uj vs. driving 
voltage, plotted for different values of 8 with v = (solid), 
and different values of v with <5 = (dashed). 



At this frequency the resistive and capacitive components 
contribute equally much to the overall cell impedance. 
Figure [7f b) clearly shows that as the voltage is increased, 
the double-layer capacitance grows, and the characteris- 
tic frequency shifts down. 

The voltage dependence of uj is shown in more detail 
in Fig. [31 where oj is plotted versus V for different values 
of 5 and v. The GCS model simply predicts lj should 
drop from uj — 1 + S at low voltage to oj ~ 5 at higher 
voltage. The same trend is seen for the Bikerman model, 
up to the point where steric exclusion sets in; beyond this 
the double-layer capacitance decreases and uj increases, 
scaling as lj = 0(V vV) at large voltage. These qualita- 
tive features predicted by our analysis may be interesting 
to compare to experimental impedance measurements at 
large ac voltages, below the threshold for Faradaic reac- 
tions or specific adsorption of ions, to seek evidence of 
steric effects in the liquid phase. 



D. Neutral salt adsorption 

In response to the ac driving, the diffuse layer period- 
ically absorbs and expels an excess amount of ions. At 
low voltage the charging comes about from both uptake 
of counterions and expulsion of coions, so the net salt 
adsorption is low, w s» q 2 /4, cf. Eq. (|56|) for q -C 1. 
At higher voltage there are essentially no more coions to 
expel, so the charging process is dominated by uptake of 
counterions and w « \q\. 

The excess amount of (counter) ions is taken up from 
both the adjacent diffusion layer and from that at the 
opposite electrode. In order to estimate when this effect 
starts to significantly perturb the concentration in the 
diffusion layer, it is necessary to know the time-average 
salt uptake (w). This is shown in Fig. Has a function 
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FIG. 9: (Color online) Contour plot of (w), the time average 
excess salt concentration in the diffuse layer, as a function of 
driving frequency and voltage for 5 = and v = 0.01. The 
dashed line marks the characteristic frequency ui . 

of driving voltage and frequency for the Bikerman model 
with 8 = and v = 0.01. At low frequency, w < w„, 
the double layer is almost fully charged, so that * « 
V cxt - At low voltage, V < 1, the figure shows that (w) ~ 
V 2 /8, while at high voltage, V > 30, the steric effects 
dominate and (w) w yjvjv. For comparison, the GCS 
model predicts (w) « 2V/ir5 in this limit, and the GC 
model (w) w exp(F/2). At high frequency, ui ^> uj , the 
bulk resistance dominates the cell impedance, so J k 
V cx t and q = 0(J/lu), from which the scaling is (w) 
V 2 /8uj 2 < 1 or (w) « 2V/-K0J > 10, depending on the 
level of charging. 

E. The limit of ionic liquids 

The weakly nonlinear regime in a blocking electrolytic 
cell generally breaks down at large voltages, when the 
neutral salt adsorption is by the diffuse layers is strong 
enough to significantly deplete the quasi-neutral bulk so- 
lution in the diffusion layers [l|, [52|, [H, This phe- 
nomenon, however, relies on the availability of available 
space in the liquid (free of ions) for the total density of 
ions to become much more concentrated in one region 
(the double layers) at the expense of another region (the 
bulk diffusion layers), which is controlled in our MPNP 
model by the parameter v = 2a?c* = 2c*/c max . In liq- 
uid electrolytes, v represents the bulk volume fraction 
of (all) solvated ions, which is typically much less than 
one, and even in saturated solutions of highly soluble ions 
would rarely exceed 0.1. As such, strongly nonlinear ef- 
fects must generally be considered (below) in electrolytes 
at large applied voltages, especially in small systems. 

The situation is different in ionic liquids or molten 
salts, which may be described by the limit v — > 1 in 
our MPNP model. This corresponds to the mean-field 



theory proposed by Kornyshev [70j where a value v < 1 
could model a somewhat lower volume fraction of the 
quasi-neutral bulk liquid phase, compared to the charged 
double layers, where strong normal electric fields may 
compress the counterions against a charged surface (as 
described above). In a molten salt, this density varia- 
tion may be comparable to the expansion upon melting 
of an ionic crystal, which can be as large as 20%, so 
we might expect v to be as small as 0.8, which is still 
much larger than for a typical electrolyte. This simple 
approach has had some success in describin g exp eriments 
and simulations of simple ionic liquids [7ll \12\ , in what 
we would call a weakly nonlinear approximation, where 
the voltage-dependent quasi-equilibrium double-layer ca- 
pacitance is coupled to a constant bulk resistor. 

An important prediction of our analysis is that this 
picture always remains valid up to large applied voltages 
for sufficiently large ^, so that ionic liquids can generally 
be described by the simple, weakly nonlinear approxi- 
mation. For a highly concentrated electrolyte (still with 
v <C 1), we can use the estimate of Kilic et al. [54| for the 
critical voltage V c (defined by e(w) = 1) to significantly 
deplete the steady-state bulk salt concentration, 

v 2(ze) 2 L 2 a 3 c* 2 ( c* V 

which grows with concentration like v (for fixed ion 
size). The basic picture is sketched in Fig.QIJl where this 
transition is represented by the dotted line. Before this 
transition is reached, the weakly nonlinear approxima- 
tion breaks down due to significant concentration varia- 
tions, which are not large enough to deplete the bulk and 
remain confined to the diffusion layers. We estimate and 
discusse these transitions below in terms of voltage, but 
here we note that they also rise steeply with concentra- 
tion in the limit of ionic liquids, v — > 1. In a molten salt 
v 1, the strongly nonlinear regime disappears, and the 
nonlinear RC circuit approximation holds for all voltages. 

The weakly nonlinear dynamics of ionic liquids in the 
MPNP model are not very different from those of con- 
centrated electrolytes at large enough voltages to trigger 
steric effects in the double layer. In Figs. O and [5] we have 
also included curves for v = 0.1 and v = 1.0. 



V. STRONGLY NONLINEAR REGIME 

The "strongly nonlinear" regime defined in Ref. [l| is 
characterized by significant 0(1) perturbations to the 
salt concentration in the quasi-neutral diffusion layers. 
The perturbations are driven by the uptake of an ex- 
cess amount of salt w = ew into the diffuse layer from 
a diffusion zone of width y^e/uj. This induces a lo- 
cal 0(ew I \Je/oj) drop in the concentration, and there- 
fore we expect the strongly nonlinear regime to start at 
V^w) = 0(1). For e = 0.001 and to = 0(1), Fig. M 
indicates this is reached for V ~ 30. 
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FIG. 10: Sketch of the different dynamical regimes for a block- 
ing cell in the space of applied voltage, V (scaled to kT/e), 
and nominal bulk volume fraction of ions, v = c*/c max . Lin- 
ear response (below the dashed line) holds for V <C 1 for any 
v and diffuse-layer thickness e = A_d / L. For thin diffuse layers 
(e <C 1) in electrolytes (y <C 1) there is a transition for V > 1 
to weakly nonlinear dynamics, where the diffuse layer acts as 
a voltage-dependent capacitor in series with a constant bulk 
resistance; at larger voltages, there is a transition to strongly 
nonlinear dynamics, which occurs first only with the oscillat- 
ing diffusion layers (dash-dot line); at higher voltages there 
is another transition (dotted line) where the bulk solution 
becomes uniformly depleted by time-averaged mass transfer 
into the diffuse layers. The transition curves rise steeply with 
v. For ionic liquids and the molten salt limit, v w 1, only 
the weakly nonlinear regime is possible, since there is not 
enough volume available to compress significant numbers of 
ions in the diffuse layers, which approach the molecular scale 
a < Ad. 



From a physical point of view, this regime of the model 
is novel and interesting in several ways. It predicts the 
possibility of "capacitive desalination" of the bulk solu- 
tion by an ac voltage, which is a remarkable example of 
rectification by nonlinearity, since even strong ac volt- 
ages are normally assumed not to perturb the bulk so- 
lution, in the absence of Faradaic reactions. This phe- 
nomenon may have interesting applications in microflu- 
idics, since ac voltages are often used to apply large elec- 
tric fields without triggering reactions. Second, concen- 
tration gradients in the oscillating diffusion layers can 
be come large enough to cause nearly complete deple- 
tion of salt just outside the double layer, causing to lose 
its quasi-equilibrium structure. This situation of "tran- 
sient limiting current" is analyzed in the next section, 
but first we describe strongly nonlinear dynamics with- 
out diffusion limitation. 
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FIG. 11: (Color online) Strongly nonlinear response at V = 
30, lu = 0.3, 8 = 0.3, v = 0, and_e = 0.001. (a) Distribution of 
the cell voltage, 14 xt = J/c a + £ — qS (dotted), onto the bulk 
electrolyte, J/c (circles), diffuse layer, £ (squares), and com- 
pact layer, — qS (triangles), (b) Concentration c s at the inner 
edge of the diffusion layer. Symbols show results from our 
full numerical solution of the PNP equations, while the solid 
lines are predictions of the (much simpler) uniformly valid 
asymptotic approximations, which are seen to be in excellent 
agreement. 



A. Dynamical response 



Figure [TT] shows the strongly nonlinear dynamic re- 
sponse at at V = 30, u — 0.3, S = 0.3, v = 0, and 
e = 0.001. First off we note that the qualitative differ- 
ence against the weakly nonlinear solution from Fig. [Sfb) 
is fairly small, even though the surface concentration c s 
shows a significant variation. Quantitatively the largest 
difference is on the zeta potential, reaching 12% rela- 
tive difference between the weakly and strongly nonlinear 
models. Perhaps this should not be too surprising: The 
surface concentration affects the double-layer charging 
dynamics only through the diffuse-layer charge-voltage 
relation, Eq. (|5ip . and only in a square- root dependence. 
Moreover, at this voltage the diffuse-layer capacitance is 
large enough that the compact layer dominates the over- 
all response. Hence for smaller values of S, we should 
see a more significant difference between the weakly and 
strongly nonlinear regimes. On the other hand, for v ^ 
the double-layer voltage eventally becomes dominated by 
the condensed phase of ions developing at the steric limit, 
which scales as |C| ~ vq 2 /2 independent of c s . 
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B. Numerical validation 

In order to test our uniformly valid asymptotic approx- 
imations above in the strongly nonlinear regime, we com- 
pare the results to the full numerical solution of the PNP 
model from Fig. [2] In Figure QT] the solid lines show the 
results from the asymptotic analysis, and symbols show 
corresponding output from the full PNP model, deter- 
mined in the following way: The compact-layer voltage, 
— <7 PNP 5, is given directly by Eq. (f2"Tj) , the bulk current, 
J PNP , and salt concentration, c PNP , are evaluated at the 
center of the cell at x — 0, the diffuse-layer voltage, i^ PNP t 
is computed as the potential drop from the electrode sur- 
face at y — (i.e., x = —1) to a point immediately out- 
side the diffuse layer, chosen (arbitrarily) at y = 3e, and 
likewise the concentration c PNP is evaluated at y = 3e. 

Overall, the agreement between the full PNP numeri- 
cal solutions and the uniformly valid asymptotic approx- 
imations is excellent, in spite of the dramatic mathe- 
matical simplification at large voltages. The bulk cur- 
rent in the asymptotic model is slightly too small when 
c s is maximal, and slightly too large when c s is mini- 
mal, with a maximal relative error of 1%, measured as 
maxt | J— J PNP |/ maxt | ^yT'NF' | _ This small discrepancy is 
primarily due to our neglect of the change in conductiv- 
ity in the diffusion layer and the associated (small) excess 
voltage, cf. Eq. f|48[) . The compact-layer voltage agrees 
very well with the full numerical solution, whereas the 
diffuse-layer voltage £ appears to be about 6% too large. 
However, the excess potential -0 m the diffuse layer falls 
off exponentially at large y, cf. Eq. (I32|) . so measuring 
the diffuse-layer voltage in the PNP model from y = 
to y = 3 we miss a (small) fraction of the "true" result. 
Accounting for this, we find the relative error is only 2%, 
mainly due to a phase lag between the two solutions. 
The same arguments apply to the salt concentration c s 
at the inner "edge" of the diffusion layer: Using Eq. (|3"0|) 
to compute c at y = 3 the agreement with the full nu- 
merical solution is accurate to within 1%, against 4% for 
the "raw" c s data in Fig. \TT\b). 

C. Local salt depletion 

In order to quantify the strength of the nonlinear re- 
sponse, we measure the minimal concentration in the dif- 
fusion layer over one period in time. For example, in 
Fig. [TlT b) the minimal concentration is about min t c s ~ 
0.5, which is attained just after t = 5.24 and again after 
t = 15.71. Figure [T2l shows the result for mint c s as a 
function of driving frequency and voltage. Figure I12f a) 
shows the result for the GCS model with 5 = 0.3, v = 0, 
and e = 0.001, with at least two important points to note: 
First, at a given driving voltage, the salt depletion is 
most significant just around the characteristic frequency 
u) , and second, at a given driving frequency min t c s falls 
off roughly linearly with V, scaling as 1 — 0(V y/ew / 8) 
for lo < L0 o . 
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FIG. 12: (Color online) Contour plot of the minimal salt con- 
centration mint Cs in the diffusion layer as a function of driving 
voltage and frequency for (a) 5 = 0.3, v — 0, and e = 0.001; 
(b) S = 0, v = 0.01, and e = 0.001. The dashed line marks 
the characteristic frequency uj and the shaded area marks 
a regime where mint c B drops to zero and the double layer is 
driven out of quasi-equilibrium. (For the lowest frequencies in 
the figure, the 0(-\/e/u) ~ 0.3 diffusion layers extend across 
most of the bulk and no longer act as mathematical boundary 
layers.) 



The frequency dependence can be understood as fol- 
lows. Earlier we argued that when the diffuse layer ab- 
sorbs neutral salt from a diffusion layer of width y^ejuj, 
the variations in c s should scale as y / etj(u>), which ex- 
plains why the salt depletion becomes less significant at 
low frequency. On the other hand, the double layer only 
gets fully charged when the system is driven below the 
characteristic frequency uj , cf. Fig. O so that overall we 
should indeed expect to see the strongest salt depletion 

for LO K L0 o . 

Another important feature of the strongly nonlinear 
regime is the possibility of transient diffusion limitation. 
This occurs when the voltage is sufficiently large, and 
the frequency sufficiently small, to temporarily, but com- 
pletely, deplete the salt concentration at the inner edge 
of the diffusion layer. The shaded area in Fig. [T2Ta) 
marks the parameter range where min t c s hits zero, and 
the quasi-equilibrium structure of the double layer breaks 
down. In this novel regime, we must revise our asymp- 
totic analysis to produce uniformly valid approximations 
accounting for transient space charge formation. This is 
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the subject of Sec. IVII below. 

Volume constraints can have a significant effect on the 
strongly nonlinear dynamics of our model problem. Fig- 
ure [T27 b) shows the corresponding results for Bikerman's 
model with 6 = 0, v = 0.01, and e = 0.001. Again, 
the salt depletion is strongest when the system is driven 
around the characteristic frequency, although this has a 
different dependence on voltage, as noted above. Fur- 
ther it is clear that when steric exclusion sets in and 
the diffuse-layer capacitance decreases, the salt deple- 
tion becomes much less significant, especially at low fre- 
quency. This effect was noted by Kilic et al [H, 
for the response to a sudden dc voltage, but its influ- 
ence on strongly nonlinear ac response is more compli- 
cated. Steric effects make the shaded area of transient 
diffusion limitation span a narrower range of frequencies, 
compared to the GCS model. However, in both models, 
the shaded area starts at roughly the same voltage for 
the characteristic frequency. 



VI. BREAKDOWN OF QUASI-EQUILIBRIUM 
DOUBLE-LAYER STRUCTURE 

As we have seen in Figure I12[ when the driving volt- 
age is increased, the salt depletion in the diffusion layer 
becomes more and more pronounced, and at some point 
the minimal concentration can even drop to zero (within 
the shaded area). At that point, the quasi-equilibrium 
structure of the double layer breaks down: The chem- 
ical potential diverges, and the effective width of the 
diffuse layer grows like 0(l/y/c^). Likewise, the quasi- 
electroneutral solution in the diffusion layer breaks down 
when the concentration approaches zero: The leading or- 
der charge density in the diffusion layer can be evaluated 
from the Poisson equation, cf. Eq. ([39)1 , 

~ _ 3/2 JdyC 
<r 

At large voltage the flux into the double layer is domi- 
nated by uptake of counterions since there are no more 
coions to expel, so that at the inner edge we have 
| J | dyc/yfe and \p\ « e 2 J 2 /c 2 . Quasi-electroneutrality 
in the diffusion layer remains a good approximation only 
as long as c> \p\ or 

c s » |eJ| 2/3 . (70) 

The breakdown of electroneutrality and concommitant 
expansion of the double layer into a non-equilibrium 
structure due to transient diffusion limitation, in the ab- 
sence of any normal flux of ions at the electrodes, is a 
novel prediction of our model which we analyze in detail 
in this section. 



A. Nonequilibrium double layer 

The breakdown of quasi-equilibrium in the double 
layer and of quasi-electroneutrality in the bulk region 
is well known for electrochemical cells driven at a dc 
Faradaic current approaching the classical "limiting" cur- 
rent [45|, [48|. At the limiting dc current, the double 
layer acquires a steady non-equilibrium structure and 
expands in dimensionless width from O(e) to 0(e 2 / 3 ), 
as first described by Smyrl and Newman [4(|. Rubin- 
stein and Shtilman later showed that a "space charge" 
region completely depleted of coions can develop at an 
electrode or ion exchange membrane when driven above 
the diffusion- limited current [471 ] . In this regime one can 
identify three sublayers within the nonequilibrium double 
layer, namely HJ] 

• An inner quasi-equilibrium layer of width 0(e) at 
the electrode surface. 

• An extended "space-charge" layer of width y a > 
0(e 2 / 3 ) that is completely depleted of coions. 

• A "Smyrl-Newman" transition layer of width 
0(e 2 / 3 ) around y = y a connecting the space-charge 
layer to the quasi-electroneutral diffusion layer. 

It is exactly the same nested boundary-layer structure 
that we see here develop in a cell driven dynamically at 
very large voltage even though the electrodes are block- 
ing with no reactions taking place, and thus no normal 
flux of ions into or out of the cell. Instead, the dou- 
ble layer is driven out of equilibrium purely by nonlinear 
electrochemical relaxation within the cell, as counterions 
are absorbed into the double layer so quickly and in such 
large numbers that bulk diffusion becomes temporarily 
rate limiting, within each ac period. We now develop 
uniformly valid asymptotic approximations for this novel 
regime. 

1. Space-charge layer 

When a negative voltage is applied on the left elec- 
trode, the space-charge layer developing is one completely 
depleted of anions, c_ = (denoting variables by a 
breve accent), while the cation concentration is nonzero, 
c + > 0, but small. The ion transport is completely dom- 
inated by migration, and the flux is determined by the 
current fed into the boundary layers from the bulk 

l -c + dJ = \J{t)\. (71) 

Substituting the Poisson equation — e 2 <9 2 = and 
integrating, we get the leading order field in the space- 
charge layer 

d v 4>=-j2\J\(y -y), (72) 
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where the integration constant y (t) is born positive and 
equal to the width of the space-charge layer. The (small) 
charge density due to the counterions in the cationic 
space-charge layer is found by differentiation, 



P : 



£+ 
2 




(73) 



and the leading order potential drop across the layer by 
integration 



$ = #0)-#j/ o ) = -f \/2\J\y 3 J 2 
3e v 



(74) 



The analysis of the opposite case, where a positive volt- 
age is applied on the electrode and a space-charge layer 
completely depleted of cations develops, is fully similar. 



2. Inner diffuse layer 

Within an O(e) distance from the electrode surface, the 
counterions remain in quasi-equilibrium with a constant 
electrochemical potential at leading order [j5Cj, i.e., for a 
space-charge layer completely depleted of anions 



M+ = <t> + log 2+ - log(l - vc+/2) 

from which the cation distribution is 

1 

c+ = -Th- 



eorist. 



(75) 



(76) 



Substituting into the Poisson equation and integrating, 
we find 



dy(, 



log 



1 



(77) 



where the integration constant k is fixed as k = 
y/2\J\y = ed y (j)(0) to match the field in the space-charge 
layer. The solution for the potential can be expressed in 
integral form as 



V 



My) 

0(0) 



log 



1 



2 C 



(78) 



The difficulty is, however, that we cannot determine the 
chemical potential /}+ by matching with the space-charge 
layer because 4> ~~ * 00 an d c+ — * for y — > oo. 

In the dilute limit c+ <C 2jv the solution can be ex- 
pressed in closed form [48l.l5ll|. Rewriting in terms of the 
excess potential V> = — <f> we obtain 



4> = 2 log [1 - (1 - e c/2 )e 



(79) 



where £ is determined by the total charge q = dy(j)(0) = 
dytp(0) + k accumulated in the non-equilibrium double 
layer, 



C= -2 log 



Q_ 

2k 



(80) 



Substituting Eq. ([79]) and c+ 
find 



-29?^ into Eq. (75]) we 



/i+ = 0(0) + 2 log(2 K ) + log (1 - e 



= (£(0) + 21og(2K)+log 



q + k 



(81) 
(82) 



We note that for q ^> k the chemical potential approaches 
a level fi+ — 0(0) + 2 log(2K) that is independent of C or 
q and determined only by the matching field from the 
space-charge layer. 

Now, provided k is much smaller than the field at the 
onset of steric exclusion, i.e., k -C \f2jv, we can use 
Eq. (JHU) for the chemical potential also for the MPB 
problem. Substituting into Eq. ([77)1 we then get 



q = ^k 2 + H log [l + 2i/k 2 (1 - e C/2) e -C] , 
and, finally, solving for <^ and using k -C \j2jv 

c 



-2 log 



+ i- v /2(^ 2/2 - 1 )A 



(83) 



(84) 



3. Transition layer 

The solution in the space-charge layer cannot be 
matched directly to the quasi-electroneutral diffusion 
layer: The leading order field from the space-charge layer 
is 0(e~ 1 ) but vanishes at y = y a , whereas the field in the 
diffusion layer is 0(1) but diverges like — J/c w l/(y~y ), 
cf. Eqs. (|45 | and (T J) . 

In the transition zone for \y — y \ < Q(e 2 ^ 3 ) the field 
has a unique profile that can be expressed as [5C| 



JP^) 

|eJ| 2 / 3 



(85) 



where P(z) is a Painleve transcendent of the second kind 
and z is a rescaled spatial variable given by 



\J\(y -_Vo) + Cs 
|eJ| 2 / 3 



(86) 



See the Appendix for more details and a plot of P(z). 

The voltage on this narrow layer is negligible in the 
overall cell response, but the solution does allow us to 
understand better both the spatial transition, and the 
transition in time from quasi- to nonequilibrium struc- 
ture. 



4- Transition from quasi- to nonequilibrium 

One challenging aspect in modelling the dynamics of 
the system is that (unlike the steady-state dc case ana- 
lyzed in all prior work) it is not sufficient to have valid so- 
lutions for quasi-equilibrium with c s S> e 2 / 3 and nonequi- 
librium with y a 3> e 2 / 3 . The dynamical response passes 
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back and forth from quasi-equilibrium to nonequilibrium, 
but the previous analysis of the inner diffuse layer from 
Sees. IIII CI and IVI A 21 breaks down in the transition 
regime from c s = 0(e 2 / 3 ) to y a — 0(e 2 / 3 ), leaving the 
diffuse-layer voltage from both Eqs. ([53]) and (|84[) diver- 
gent for c s —>■ and y — ► 0, respectively, at fixed q. This 
is problematic since the charge-voltage relation plays a 
central role in our dynamical model. 

In order to resolve this, we have developed an approx- 
imate solution of the standard PNP equations for the 
inner diffuse layer, that is uniformly valid both in quasi- 
equilibrium, nonequilibrium, and across the transition 
regime. Essentially, our approximation amounts to as- 
suming constant, rather than variable, coefficients in the 
equation for the excess field in the boundary layer from 
the Painleve II problem; see the Appendix for technical 
details. Then we obtain 



-0 = 4 tanh 



tanh(C/4)e- K ^ 



1 - *tanh(C/4)(l- e" K S) 



(87) 



where 



K = \eJ\^ 3 ^P(-z ) 2 /2-z , 
k= |eJ| 1/3 sign(J)F(-z ), 
z = (\j\yo-c s )/\eJ\ 2l \ 

and the charge- voltage relation is 



(89) 
(90) 



C = -21og 



q + k 



2(k + k) 



K + K 2(k + K) 



(91) 



The asymptotics of P(z) are such that in the quasi- 
equilibrium limit, where c s 3> e 2 / 3 and z <C — 1, we 
have P(— z a ) w l/z a , k Vc~ s , k « 0, and we recover 
the standard Gouy-Chapman solution. In the nonequi- 
librium limit, where y 3> e 2 ^ 3 and z 3> 1, we have 
P{— z Q ) ps —y/2z , k as \k\ w y/2\J\y , and we recover 
the result from Eq. ([80]) . 

Extending this analysis to account for steric exclusion 
does not affect the space-charge layer, since the concen- 
tration there is low. It does, however, affect the solution 
in the inner diffuse layer. We are not able to produce 
an explicit analytical solution like Eq. (|87|) to the MPNP 
equations. But, noting that the MPB charge-voltage re- 
lations both quasi-equilibrium and nonequilibrium cases, 
Eqs. (f53|) and {84|) are obtained by substituting q in the 
corresponding PB results with sign(q)y/2(e l/ Q 2 / 2 _ V)Jv. 
Extrapolating this observation to the general form in 
Eq. (|9"Tj) , we argue that the general MPB charge volt- 
age relation is simply obtained by replacing q with 
siga{q)y/2(e v ^ 2 / 2 - l)/v in Eq. 



(j9lj) . At least this form 
reduces to Eq. {91]) for q 2 < 2/v, to Eqs. ([53]) and $Q 

for c s 3> e 2 / 3 and y a e 2 / 3 , respectively, and to 
£ w — sign(<7)^<7 2 for q 2 3> 2/v, 



B. Modified diffusion layer 

The width y {t) of the space-charge layer is equal to 
the width of the region of complete salt depletion in the 
diffusion layer. In the quasi-electroneutral part we still 
need to solve a simple diffusion problem 

d t c = d 2 c, 

however, now it is to be solved on the dynamically chang- 
ing interval y € [y (t] >oo), where y = y /\f^ is deter- 
mined such that 



Vo(t) 

c{y {t),t) 



for 
for 



6(0, t) > 0, 
Voit) > 0. 



(92) 
(93) 



At the inner edge of the diffusion layer, the boundary 
condition is still obtained by matching with the salt flux 
out of the double layer 



^= lim dy£ = lim F = F (t). 

V 6 y^io y^iia 



(94) 



The problem can be solved using the method of images 

c = c + V^J q \[Gu(\y~y (t%t-t') 

+ G u (y + y (t'), t - t')]F (t') dt', (95) 

where the injection of flux at both ±y (t') makes the ori- 
gin act like a reflecting boundary, and the Green's func- 
tion G^y.t) is defined in Eq. ([4*6"j). The definition of y Q 
through Eq. ([93)) ensures that 6 = and d y 6 — for 
V < Vo{t), and hence, using Td y G LJ (0 ± ,t) = S + (t), the 
boundary condition ([9"4"| is indeed satisfied. 

Interestingly, the modified diffusion layer does have 
a time-average excess salt content relative to the bulk: 
Substituting Eq. ([H in Eq. {95]) we find that 



(c) =c -\((\y- Vol + \y + Vo\)F ), (96) 

and hence the time-average excess amount of salt con- 
tained is 

f°° e 
V~e(w) = V~e (c)-cdy = --(y F ), (97) 



where we used that (F Q ) = -(d t w) = [l04j . Note that 
since F is generally negative when y D is nonzero, the 
excess salt content will be positive. 

Since the local conductivity goes to zero when 6 van- 
ishes at the inner edge, one might worry that the poten- 
tial drop across the diffusion layer could be large. How- 
ever, because the quasi-electroneutral solution breaks 
down for c < |eJ| 2 ^ 3 , we find that the overall ohmic 
potential drop across the quasi-electroneutral diffusion 
layer is limited to 0( log |e J| 2 / 3 ) , which remains a small 
perturbation to that across the bulk electrolyte. 
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FIG. 13: (Color online) Strongly nonlinear response at V = 
120, w = 0.3, S = 0.3, v = 0, and e = 0.001. (a) Distribution 
of the cell voltage, V ex t = J/c + & + ( — qS (dotted) , onto the 
bulk electrolyte, J/c (circles), diffuse layer (dashed), space- 
charge layer (dash-dot), and compact layer, —q~8 (triangles); 
also shown is the sum of the space-charge and diffuse-layer 
voltages, £ + $ (squares), (b) Concentration c s at the in- 
ner edge of the diffusion layer (triangles) and extent y a of 
space-charge layer (crosses). Solid and broken lines show re- 
sults from our asymptotic model, and symbols show the full 
numerical solution of the PNP equations. 



C. Model summary 

Summarizing the model for the leading order dynamic 
out-of-quasi-equilibrium response, the charging of the 
double-layer at the leftmost electrode is governed by 



V cxt = J/c + <i> + (-qS, 
d t q = — J, 

w = \Jq 2 + 4c s - y/4£a, 
d t w = -F Q , 



(98) 
(99) 
(100) 
(101) 



with the diffuse- and space-charge layer voltages Q and <& 
given by 



C=-21og 



Q + K 



2(k + k) 



Q + k 



K + k 2(k + k) 



$ = sign(J)f sj2\J\yl' 2 
3e v 



(102) 
(103) 



Here Q = sign(q)y/2[exp(vq 2 /2) - \\jv in the MPNP 
model, reducing to Q = q in the standard PNP model. 



The salt concentration c s at the inner edge of the diffu- 
sion layer is determined from 



. 1 



c s = c + ^/e— 



1 



G u (\y (t)-y°(t')\,t-t') 
+ G w {y (t) + y (t'),t - 0] F {t r ) di', (104) 

and the width y Q = y/ey of the space-charge layer from 
c s y o = 0, c s >0, y o >0. (105) 
The parameters k and k are given by 

k = \eJ\V 3 y/3P(-z )*/2-z , (106) 
K = sign(J)|eJ| 2 / 3 P(-z ), (107) 



where P{z) is the Painlcve transcendent and z Q = 
{\J\Vo — c s )/|eJ| 2 / 3 is the rescaled position of the transi- 
tion layer relative to the electrode. Finally, the bulk con- 
centration is again determined by imposing global con- 
servation of salt in the cell to get 



c = l + ^(y F )/2-e(w). 



(108) 



We solve this problem numerically using a timestep- 
ping algorithm. The major difficulty is to determine y (t) 
in a self-consistent way, which we achieve using a bisec- 
tion algorithm. More details are given in Ref. [9fJ. 



D. Dynamical response 

Figure [T^] shows the strongly nonlinear dynamic re- 
sponse at V = 120, u) — 0.3, S = 0.3, v = 0, and 
e = 0.001. In particular we note in Fig. [TSTb) the ap- 
pearance of a transient space-charge layer extending to 
a width of y < 0.04, with an associated voltage <& in 
Fig. [TSTa) that induces a visible drop in the bulk current 
J. 

In order to validate to asymptotic analysis we compare 
the results to the full numerical solution from Fig. [3l 
as shown with symbols in Fig. [T31 Like in Sec. IV Bl 
the compact layer voltage — g PNP i5 is given directly by 
Eq. (|21j). and the bulk current J™ 13 and salt concentra- 
tion c PNP are evaluated at the center of the cell. The 
width y PNP of the space-charge layer is taken (arbitrar- 
ily) as the largest region where either the cation or anion 
concentration drops below i|eJ| 2 / 3 , and, finally, the con- 
centration c PNP and the overall voltage (C + $) PNP across 
the diffuse and space-charge layers are evaluated at the 



position ?/ PNP + 3e. 

The agreement between the asymptotic approximation 
and the full numerical solution is rather good, although 
the bulk current J is visibly somewhat too low (by w 5%) 
when c s is large and too high when the c s is small. As 
before in Fig.[TTl part of the discrepancy on l>, and c s 
is due to the difficulty with defining the diffusion-layer 
inner "edge" on the full numerical solution, so it may be 
more appropriate to compare the full spatial profiles. We 
proceed to do that in the following section. 
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FIG. 14: (Color online) Concentration profiles at V = 120, ui = 0.3, S = 0.3, v — 0, and e = 0.001; "frames" cover one-half 
period in time starting from t = 0. Open and filled circles show the cation and anion concentration, respectively, according to 
the full numerical solution from Fig. [3] and solid lines show uniformly valid approximations based on the asymptotic analysis, 
with very good qualitative and quantitative agreement. Black arrows mark the extent y (t) of the space-charge layer (if any) 
according to the asymptotic model. 
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FIG. 15: (Color online) Strongly nonlinear response for MPB model with relatively large and small bulk ion volume fraction, 
v = 0.01 and 0.0001, respectively, at V = 120, uo = 0.3, S — 0, and e — 0.001. (a) and (c) Distribution of the cell voltage, 
Vcxt = J/ c + & + C (dotted), onto the bulk electrolyte, J/c a (solid), diffuse layer (dashed), and space-charge layer (dash-dot), 
(b) and (d) Concentration c 3 at the inner edge of the diffusion layer (solid) and extent y of space-charge layer (dashed). 



E. Uniformly valid approximations 

So far we have been focusing on integral quantities 
such as the total charge and voltage, but the asymp- 
totic analysis also predict the full spatial profiles for the 
potential and ion concentrations in the cell: Uniformly 
valid approximations in space are constructed by adding 
the inner and outer approximations and subtracting the 
overlaps [J, [54| . In the absense of a space-charge layer 
the ion concentrations are given by 



c± (x,t) 



l + x 

C± I , t I - C, 



l + x 
c I — =-J 



1 — x 



(109) 



In the presense of a space-charge layer, the challenge is 
to tie up the O(e) but divergent countcrion concentration 
in the space-charge layer, cf. Eq. ([75)1 . with the 0(1) 
but vanishing concentration in the diffusion layer at y = 
y . The key is to employ the solution from the Smyrl- 
Newman transition layer: In the transition to the space- 
charge layer on the left electrode, the concentration can 
be written as 



c± = \eJ\ 2/3 (z + P 2 /1 T sign( J)d z P), 



(110) 



as discussed in the Appendix. For y ^> y Q this reduces 
to c + = c_ = {y — y )\J\, matching the flux d^cj\fi w 
|J| at the inner edge of the diffusion layer, whereas for 
y <C y a the space-charge density Eq. ([75)) is recovered. It 



is convenient to rewrite Eq. (|110[) as c± = c + <f ± where 
<f± = c± — c is the excess concentration in the space- 
charge and transition layers relative to the diffusion layer, 
given by 

£t = |eJ| 2 / 3 (min{z, 0} + P 2 /2 T sign( J)d z P). (Ill) 

In line with this, the excess concentration in the inner 
diffuse layer can be expressed as 

ftfc = \eJ\ 2/3 (R 2 /2 + PR =F sign(J)d z R) 

= (dyi>) 2 /2 + ndy4> T dli> (112) 

where R = sign(J)d z rp = sign( J)dy^)/\eJ^' 3 is the 
(rescaled) inner excess field. Substituting Eq. ([57)) we 
obtain the following lenghty expression 



=F 4k 



2 2S - k sinh(C/4)(ee K S - Se^ /6( 
[2Ksinh(C/4) + Qe K v - Ee~ K y /6] 2 
2 Ksinh(C/4)(ee K i' + Se-^/e) 
[2Ksinh(C/4) + 6e re « - Se-^/G] 2 ' 



(113) 



where S and are shorthands for 5 = (k 2 — 
K 2 )sinh 2 (C/4) and 6 = «:cosh(C/4) - Ksinh(C/4). With 
this, the general form of (our approximation for) the 
ion distributions, uniformly valid in space, and in time 
from quasi-equilibrium, across the transition regime, to 
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nonequilibrium, becomes 

c±(x,t) =?± +Z±{z+(t),t) 

+ « F («_(t),t)+ft F ^i^£,t^ ) (114) 

where z±(i) = [|J|(1 ± x T Vo{t)) ± c a (t)] /|e J| 2 / 3 . 

The resulting concentration profiles are shown in 
Fig. [H for the solution at V = 120, u = 0.3, S = 0.3, 
2/ = 0, and e = 0.001, displaying first injection of salt 
from the double layer into the diffusion layer, followed by 
re-uptake un the double layer, salt depletion in the dif- 
fusion layer with formation and growth of an extended 
space-charge region, and, finally, collapse of the space- 
charge layer when the cell current changes direction. The 
figure also compares the uniformly valid approximation 
to the full numerical solution from Fig. [3l and shows 
very good agreement. The relative error, measured as 
\c±(y,t) - c^ NP (y,t)\/max t c^ NP (y,t), is below 4% for 
all y on the eight frames displayed in Fig. [T4l 

F. Dominant voltage in double layer 

For electrochemical cells running with a dc Faradaic 
current, it is well known that concentration polarization 
can play a dominant role at large voltage, with the space- 
charge layer determining the overall current- voltage re- 
lation for the system [43, UK]. It is clear from Fig. [T3] 
that although $ is smaller than the compact-layer volt- 
age, —qS, it does affect the bulk current and slows down 
the charging process. 

On the other hand, since $ depends on the capaci- 
tive current in ac, and since the diffuse-layer capacitance 
eventually drops when crowding starts to kick in, one 
might ask if the overall cell response will not be domi- 
nated by £ at large voltage? Of course, that will depend 
on just how early the steric limit is reached, i.e., it de- 
pends on the nominal ion volume fraction v. 

Figure [15] shows the strongly nonlinear response for 
the Bikerman model with two different values of v, at 
V = 120, u = 0.3, 5 = 0, and e = 0.001. In Fig. E^a) 
and (b) where v = 0.01, the diffuse-layer voltage £ dom- 
inates in the cell while the space-charge layer voltage <J> 
is negligible. Also the bulk ohmic potential drop J/c a 
is small because the driving frequency is below the char- 
acteristic frequency, w < uj , cf. Fig. [5] In Fig. rTSTc) 
and (d) where v — 0.0001, the situation is the opposite: 
$ dominates over £ in the non-equilibrium double layer, 
while J I c dominates the overall cell response because 
this system is driven above the characteristic frequency, 

u> > L0 o . 

The competition between $, £, and — qd is investi- 
gated further in Fig. 1161 showing the peak values max t 



maxj (, and max t q5 as a function of driving voltage V 
in panels (a) to (c) and frequency u> in panels (d) to (f). 
Figure UM&) and (d) shows results for the PB model with 
5 = 0.3. Here the compact-layer voltage dominates, al- 
though $ grows to a significant fraction at large voltage. 
Note in Fig. \W[d) that $ oc 3 1 l 2 y ? 'J' 1 peaks around the 
characteristic frequency oj « 0.3: At higher frequencies 
the double layer is not fully charged so w and y de- 
creases, while at lower frequencies it is fully charged and 
hence J decreases. 

Figure RTOf b) and (e) shows results for the MPB model 
with v = 0.01, where £ completely dominates over 
whereas in Fig. fTHT c) and (f) with v = 0.0001, $ domi- 
nates over £ at large voltage and not-so-high frequency. 

The relative magnitude of the double-layer voltages 
can be understood from a simple estimate: Once steric 
effects dominate in the diffuse layer we have 

C = 0(uf/2). (115) 

For the space-charge layer we have J = 0(ujq) and y a = 
0(ew) — 0{eq) so that 

* = 0(| J\ 1/2 y 3 /2 /e) = 0(V^q 2 ). (116) 

This is an important finding: With both £ and $ scaling 
as 0(q 2 ) at large voltage, we expect C to dominate over 
$ for v 3> \/ecJ, i.e., in systems with high nominal con- 
centration (large v), large electrode separation (small e), 
and at low frequency (ui <C lo ). Conversely, we expect $ 
to play a dominant role at large voltage for systems with 
very dilute electrolytes and small (micro) electrode geom- 
etry, driven around the characteristic (RC) frequency. 

VII. SUMMARY AND DISCUSSION 

We have developed a dynamical model for the re- 
sponse of dilute electrolytes to large applied ac voltages, 
building on a body of theoretical work on diffuse-charge 
dynamics for both the weakly and strongly nonlinear 
regimes [l], [52l IBij , and on concentration polarization and 
space-charge layers in dc electrochemical systems run- 
ning at steady-state conditions [13, [H, [5(|. Our origi- 
nal contributions are the solution in the oscillating diffu- 
sion layer, controlling the extent of the transient space- 
charge layer, and the uniformly valid formulation of the 
charge- voltage relation over the transition between quasi- 
equilibrium and non-equilibrium. 

We have compared our asymptotic analysis for the 
PNP model to a full numerical solution of the PNP equa- 
tions, and found good qualitative and quantitative agree- 
ment. The strongly nonlinear regime, characterized by 
strong concentration gradients in the diffusion layer, set 
in for y/euj(w) = 0(1), where the time- average excess salt 
concentration (w) depends on the driving frequency and 
voltage, but also on the intrinsic surface capacitance and 
crowding effects through S and is, respectively. At very 



24 




FIG. 16: (Color online) Peak voltages on space-charge layer, max; $ (solid), diffuse layer, maxt ( (dashed), and compact layer, 
maxt qS (dash-dot), for different values of the capacitance ratio 5 and nominal ion volume fraction v. Panels (a) to (c) show 
results at ui = 0.3 as a function of V, and (d) to (f) show results at V = 120 as a function of ui. For the PB model, (a) and 
(d) show that the compact layer voltage — qS dominates at S — 0.3, although 4> becomes significant at large voltage. For the 
MPB model, (b) and (e) show that the diffuse-layer voltage ( dominates at v — 0.01, while in (c) and (f) the space-charge 
layer- voltage $ dominates at v — 0.0001. 



large voltage we argue that the cell response should be 
dominated by space-charge for v <C y/ao and by crowding 
effects for v ^> y/euj. 

Recently, Beunis et al. [55| presented an analysis of 
the transient response to a dc step large enough to in- 
troduce transient space-charge. They analyse four ex- 
treme cases: The "double-layer limited" (V <C 1, e <C 1), 
"diffusion limited" (V <C 1, e ^> 1), "geometry lim- 
ited" (F > 1, e > 1/VV), and "space-charge limited" 
(V >■ 1, e <C 1/vt), and develop closed form analyt- 
ical solutions in each of those limits. In particular, for 
the space-charge limited response, setting $ = V and 
Vo = e l<zl> they predict a characteristic 0(t~ 3 ^ 4 ) depen- 
dence in the bulk current, which they verify by experi- 
ments on a system with surfactant micelle charge carriers. 



While the simple analytical results of Ref. [5|| provides 
important insight to the limiting case when the space- 
charge layer completely dominates the response, our dy- 
namical model is uniformly valid from small to very large 
voltage. The general applicability, however, comes at the 
expense that our strongly nonlinear model requires a set 
of integro-differential-algebraic equations to be solved nu- 
merically. In comparison, the weakly nonlinear "circuit" 
model, which neglects any perturbations to the bulk and 
diffusion-layer electrolyte concentration at leading order, 
can be formulated as a simple ordinary differential equa- 
tion. 

We also consider steric effects of finite-sized ions in 
blocking cells under large ac voltages. This leads to a 
novel and strong dependence on the bulk volume fraction 
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of ions is, which acts as a third dimensionless parameter, 
along with e and V, to determine different dynamical 
regimes. We consider the regime of ionic liquids v = 
O(l), up to the molten salt limit v rs 1, and find that 
strongly nonlinear regime disappears with increasing v 
due to dominant steric effects, which prevent the double 
layers from adsorbing significant numbers of ions from the 
bulk. The classical diffuse layer is effectively replaced by 
a molecular condensed layer. As a result, we justify the 
use of weakly nonlinear circuit models, as in Refs. [7(| 
l7ll W\ to describe the dynamics of ionic liquids up to 
very large, time-dependent applied voltages. 

In many cases, even in dilute electrolytes, the nonlinear 
circuit model may actually give good account for over- 
all cell current-voltage response, in particular when the 
double-layer capacitance is dominated by the compact 
layer, or when crowding effects set in and a condensed 
layer forms at the electrode. However, this does not nec- 
essarily mean that the weakly nonlinear analysis will ac- 
count well for all aspects of the electrokinetic response, 
such as ac-electroosmotic fluid motion and pumping. 

A. Two or more dimensions 

The boundary-layer analysis can be easily extended 
to higher dimensions, provided the electrode geometry 
is smooth enough to be considered locally flat on the 
boundary-layer length scale. Then, the steady-state bulk 
response becomes 

V-J = and V-F = 0, (117) 

where J = —cV<t> is the current, F = — Vc+ Pe(u)c is 
the salt flux, and c = c(r) is constant in time but not 
in space. The last term in the flux describes advection 
by the average fluid velocity (u); Pe = uqL/D is the 
Peclet number, where u = e(kT j ze) 2 j r\L is the elec- 
troosmotic (EO) velocity scale, and r] is the dynamic vis- 
cosity. In the surface conservation laws, tangential flux 
through the highly charged diffuse layer must be taken 
into account [H, HE HH lH, [H[ , leading to 

d t q = n ■ J + eV s ■ J s , 
d t w = F + eV s -F s , (118) 

where V s is the tangential gradient, and J s and F s are 
the surface excess current and salt flux, respectively, due 
to surface migration and EO convection. For the PB 
model it can be shown that 

J s = (1 + Pe){wV s 4> + qV s logc s ), 

F s = {l + Pe)(qV s 4> + wV s logc s ). (119) 

The same results also apply for the MPB model, pro- 
vided the ion mobility in the highly crowded double layer 
is equal to that in the bulk, which is, however, question- 
able nil- 



Assuming transverse convection is weak enough, 
ePe|u| <C 1, the diffusion layer can still be modelled by 
simple ID diffusion in the normal direction, with the con- 
centration given by Eq. (|9"5"j). Matching with the steady 
solution in the bulk is then obtained by 

n F = (F ) = e(V s -F s ), (120) 

since the oscillating diffusion layer does not accumulate 
any salt on time average, but only acts as a buffer zone 
for the periodic flux in and out of the diffuse layer. In 
this way, surface conduction can drive bulk concentration 
gradients even in the steady-state response [Hj]. 

The bulk fluid motion is driven primarily by EO slip 
from the boundary layers. For the quasi-equilibrium dou- 
ble layer, the effective tangential slip velocity according 
to PB theory becomes 

u s = CV s + 41og(cosh(C/4))V s logc s . (121) 

For induced-charge electroosmosis (ICEO) where both £ 
and 4> depend on the external driving voltage, the ve- 
locity scales as 0(V 2 ) at low voltage [15|,[95j]. At larger 
voltage, PB theory predicts stall of £ and scaling only 
as O(VTogV) [ll[, whereas the MPB model predicts a 
return to the 0(V 2 ) scaling once crowding effects set 
in Or, this assumes the viscosity in the highly 

crowded double layer is equal to the bulk value; if it is 
significantly reduced one might expect the velocity to 
scale as 0(Vlog where log y/2/v is the diffuse 

layer voltage in the dilute part outside the condensed 
layer (57| . 

When the double layers are driven out of quasi- 
equilibrium, they are still governed by surface conser- 
vation laws like Eq. (| 1 1 8[) . The surface fluxes J s and F s 
are no longer given by Eq. (|119j) . but they remain domi- 
nated by the inner diffuse layer, since the concentration 
(and hence conductivity) in the space-charge layer is low. 

If the space-charge layer does not contribute much to 
the surface fluxes, it plays a major role on EO fluid mo- 
tion: The voltage 5> drives a Smoluchowski-type slip ve- 
locity 

u s = $V^, (122) 

for which Dukhin and co-workers coined the term "elec- 
troosmosis of the second kind" (E02) to distinguish 
it from the quasi-equilibrium response [96j . This phe- 
nomenon has been studied extensively in the context 
of nonlinear electrophoresis of conductive particles made 
from ion-exchanger material [96l |97| . 

Since many ac electrokinetic experiments involve mi- 
cro elect r odes and applied voltages of a few volt [IE [H, 
HH, llOOl . Il0l| . includin g investi gations on ac electroos- 
motic micropumps [15. Il4 Il02| , and since our analysis 
has shown that this is enough to create strong concentra- 
tion polarization and transient space-charge around the 
electrodes, we believe that E02 could be important for 
interpreting the experimental results. 
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Rubinstein and Zaltzmann showed that E02 renders 
linearly unstable the quiescent solution of concentration 
polarization on a planar permselective membrane run- 
ning at dc, leading to spontaneous formation of vortex 
pairs that stir up the concentration profile in the diffu- 
sion layer, which in turns enables the passage of "super- 
limiting" current through the membrane pOl . l9l| . Pre- 
sumably, a similar instability could occur for transient 
space-charge layers, although the threshold voltage may 
depend on whether the space-charge layer voltage $ is 
dominating in the overall cell response or not. 

B. Diffusive dynamics in the bulk 

The transient response in the bulk, while the cell 
relaxes towards the steady-state periodic solution, or 
as arising from a slowly varying ac voltage amplitude 
Kact = V (t) sin(ujt) + V (t), is governed by diffusive dy- 
namics on the slow time scale t = et 

die = —V ■ F = V 2 c — Pe(u) • Vc, (123) 

driven by the time-average flux into the double layer 

n-F = e(V s -F s )~dt(w). (124) 

Here the time averages are taken on the RC time scale 
t, i.e., for each timestep taken on the slow time scale, we 
require the periodic response of the boundary layers on 
the RC time scale to be determined. 



C. Faradaic reactions and general electrolytes 

Much interest on electrokinetics is of course associated 
with electrochemistry and reactions on electrodes that 
are not blocking but support the passage of a Faradaic 
current. Then the surface conservation laws are enriched 
(is this the right word??) by the injection of a Faradaic 
current J ox t at the electrode surface, and an associated 
salt flux F ext . 

Depending on the charge-transfer resistance and on 
the driving frequency and voltage, the Faradaic current 
may be small compared to the capacitive current in ac, 
or it may completely dominate the charging dynamics 
of the double layer. For reactions controlled by Butlcr- 
Volmcr kinetics, where the reaction rate grows exponen- 
tially with (compact-layer) voltage, the latter should be 
the case at sufficiently large voltage. Of course, this im- 
plies that the RC time scale, formed by the bulk ohmic 
resistance and double-layer capacitance, may not be ap- 
propriate for describing the system response. 

Moreover, for many reactions it is not sufficient to as- 
sume a binary electrolyte since neutral reaction prod- 
ucts also play an important role. The general problem 
of the dynamic response for general electrolytes seems 
daunting, and is a challenge even for the steady-state 
response [l03j |. 
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APPENDIX: 

The theory for non-equilibrium double layers was orig- 
inally developed for electrochemical systems passing a 
dc current in steady state [U Eg], with boundary 
conditions representing either normal flux of ions into a 
permeable electrodialysis membrane [47| or via Faradaic 
charge-transfer reactions at an electrode If we con- 
sider a cationic space-charge layer formed with negative 
voltage on the left electrode, the Nernst-Planck equa- 
tions become 

F + = -d y c+ - c+d y (f> = 2 J < 0, (A.l) 
F_ = -d v c- + c-d y (t) = 0. (A.2) 

For a dc electrochemical system, this holds across the 
entire cell. For an ac system driven at large voltage, 
this holds in the inner diffuse, space-charge, and "Smyrl- 
Newman" transition layers because the charging process 
is dominated by uptake of counterions (cations) rather 
than expulsion of coions (anions), such that |F + | ^> \F—\. 
However, in the bulk region and in the diffusion layer we 
have primarily ohmic transport and F + w — F— ~ J. 
After some manipulations on Eqs. (| A. 1 1) . (|A.2|1 . and © 
the problem is reduced to a single master equation for 
the electric field [II] 

e 2 8 2 y E - \e 2 E? - \J\(y - y* )E = |J|, (A.3) 

where E — ~d v <f> and y* is an integration constant. 
For y* > it is equal to the width y of the space- 
charge layer, whereas for y* < we interpret it as 
y* = —c„/\J\ (ij]. Rescaling with 

ITI1/3 2/3 

e = -^jt p and y-y*o = jjp>z (a.4) 

we then arrive at 

d 2 P=^P 3 + zP+l (A.5) 

This is an instance of the second order ordinary differen- 
tial equation with Painleve property (i.e., all movable sin- 
gularities are poles) defining the Painleve transcendents 
of the second kind. The connection between steady dc 
current in electrochemical cells and Painleve transcen- 
dents was first noted by Grafov and Chernenko [4l|, [12] • 
Eq. (|A.5|) has a unique "transition layer" solution P(z) 
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FIG. 17: (Color online) Solution in "Smyrl-Newman" transition layer in terms of Painleve transcendents, (a) Rescaled electric 
field P(z) (solid) and leading order terms from asymptotic expansion in space-charge layer, z <C —1, and diffusion layer, z ^> 1, 
(dotted). The dashed line shows a solution P to Eq. (|A,5j) with b.c. P — P = — 10 applied at z = —z = —5. (b) Potential 
variation <f> — — J* P(z')dz' in transition layer (solid), and leading order terms from asymptotic approximation (dotted). Again, 
the dashed line shows the result for the solution P. (c) Rescaled charge distribution d z P (solid) and leading order terms from 
asymptotic approximation (dotted). Symbols show individual ion concentrations, cf. Eq. (|A.9|) . (d) Excess voltage ( as a 
function of z with P — —10 (solid), and leading order terms from asymptotic approximation (dotted). The circle marks 
z = 5 corresponding to P from panels (a) and (b), and the dashed and dash-dot lines show the approximation by Eqs. (|A.19p 
and (IA.24[) . respectively. 



with no poles on the real axis and the following asymp- 
totic behaviour 

P(z) = ( -U^Sn/ul f° r Z " + °°' (A-6) 

[ — v— 2z + 0(1/ z) lor z — > — oo. 

The detailed shape of P(z) is shown in Fig. flTT a) and 
compared with the leading order asymptotics. Fig. 117( b) 
shows the potential variation 



0(z) - - / P(z')dz', 
Jo 

and Fig. \T7\c) displays the rescaled charge density 
P -Q P 

and individual ion concentrations 

c± - Z +-P 2 ±d P 



(A.7) 



(A. 



(A.9) 



Boundary layer 

The form of P(z) describes the solution in the interior 
of the electrochemical cell. However, P(z) generally does 
not satisfy the boundary conditions at the electrodes con- 
fining the cell. Imposing b.c.'s on the solution gives rise 
to boundary layers, that can be understood mathemat- 
ically as originating from poles in the solution, located 
outside the domain of the physical cell. 

We focus on the behaviour at the left electrode, and 
consider a solution P(z) to Eq. (|A.5[) on the interval z £ 
[— z ,oo), with boundary conditions P(— z a ) — P and 
P(oo) = 0. Here 



\J\Vo 



|eJ| 2 / 3 



(A.10) 



corresponds to the rescaled position of the electrode. Fig- 
ure [T7fa) shows the result for z — 5 and P Q = —10: The 
excess field is rapidly screened out, and for z > —4 we 
see that P(z) follows P(z) closely. 

Let us introduce the excess field R = P — P in the 
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boundary layer. Substituting into Eq. (|A.5[) we obtain 

d 2 z R=-R i + -PR 2 + k 2 R 1 (A.ll) 
z 2 2 



where k = ^/3P 2 /2 + z. Since the boundary layer is thin, 
it is reasonable to approximate the variable coefficients 
with constants P = P(— z Q ) and k Q = \J"iP 2 /2 — z Q to 
get 



d 2 R = ii? 3 + 3 -P Q R 2 + k 2 a R. 



(A.12) 



This approximation is crudest for z close to zero, where 
the local screening length l/k Q has a maximum; for z a <C 
— 1 and z„ > 1 we have k Q f=a y/—z and fc « — P D f=a 
y/2z , respectively. 

Integrating twice on Eq. (|A.12|) we obtain 



d z R = -RJR 2 /A + P R + kl, 



(A.13) 



and 



z + z Q = — 



from which 

R : 



sinh 1 



sinh 



PoR + 2k 2 



\R\y/P*/2- Zo J 

( P R + 2k 2 



\R \^/P 2 /2-z 0/ 



(A.14) 



2kl 



where 



Pa - sign{R )^P 2 /2- z a sinh [k a (z + z*)] ' 

(A.15) 



, + fsnW P -^±M ). (A.16) 

fco \\R \^/P 2 /2-zJ 



Finally, the excess potential -0 
grating R = —d z tjj to 



is found by inte- 



•ijj = 4 tanh 1 



tanh(C/4)e 



-fc (z+z ) 



f + tanh(C/4)[f - e -fco(z+z )] 



(A.17) 



2 log 



1 + |^ e C72 + (1 _ e i/2y-k ( z +z ) 



1 + [ e C/2 + (i _ eC/2) e -*.(»+ io )] 



(A.18) 



Here £ = ip(—z ) is determined through the boundary 
condition P, = P a + R a as 



C = -21og 



2 sinh 



P + P> 



i 2 



2(fc - P ) 



k Q P Q 2(/c -fo) i 
(A.19) 



log 



(A.20) 



Using k a « \J—z Q and P G « l/z Q for z G <g; — I it is easily 
verified that Eq. (|A.19p reduces to Chapman's formula 



P, 



C w 2 sinh" 1 



in quasi-equilibrium, and, similarly, using k 
\f2z~ for z„ > 1 we find 



C~-21og 



f 



Pn 



(A.21) 

-Pa « 

(A.22) 



in accordance with Eq. (|80[) in non-equilibrium. How- 
ever, while Eqs. (|A.21|) and (|A.22j) diverge for z Q — > 
at fixed P OJ our general result (|A.19[) only displays a lo- 
cal maximum. This is shown in Fig. II7f d) where C is 
plotted as function of z for P Q = —10. The figure also 
compares our result to the excess voltage from a direct 
numerical solution for P, and it is seen that Eq. (|A. 19|) 
slightly overestimates £ for z Q < and underestimates it 
for z a > 0. 



Singular transcendents 

Zaltzman and Rubinstein [5l| systematically studied 
the transition from quasi-equilibrium to noncquilibrium 
by considering various ranges for the parameter z a , solv- 
ing appropriate approximations to the Painleve equation 
in each range. In the transition regime they approach 
the innermost part of the inner diffuse layer by an alge- 
braically decaying solution, matched to a singular solu- 
tion of the full original Painleve equation (|A.5[) in the 
outer part, writing 



P = 



z + z - 2/P z + z 



+ Pt(z;z ). (A.23) 



Here 2/ (z + z Q ) + pt(z; z a ) is the regular part of a solution 
pt(z;z D ) with a simple pole at —z a , P^(z;z Q ) ~ — 2/(z + 
z ) for z — * — z G , and Pt(oo;z ) = 0. This allows them 
to write 



C = -21og(-P ) + <^(z ) 



(A.24) 



where 

ip(z ) = lim 



P t (z';z ) + - P(z')dz' 



log(z + z ) 



(A.25) 



Their approximation is highly accurate for z Q in the tran- 
sition range and large enough P a -C — 1, and it is simple 
to evaluate once the function ip(z ) has been tabulated. 
For comparison, our result (|A. 19|) has a finite error for 
z Q in the transition range, an error that does not vanish 
at large P but tends to a finite value, essentially being 
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due to our approximation of the variable screening "con- quasi-cquilibrium and noncquilibrium limits, allowing us 
stant" k(z) in the tail of the excess field by a real constant to use a single formulation of the charge- voltage relation 
k Q . On the other hand, our result matches fully with the for the entire dynamic solution procedure. 
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